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The  parabolized  stability  equations  (PSE)  axe  a  new  and  more  reliable  approach 
to  analyzing  the  stability  of  streamwise  varying  flows  such  as  boundciry  layers.  This 
approach  has  been  previously  validated  for  idealized  incompressible  flows.  Here,  the 
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to  permit  the  anadysis  of  high-speed  boundary-layer  flows  over  fairly  general  bodies. 
Vigorous  numerical  studies  are  carried  out  to  study  convergence  and  accuracy  of  the 
lineaj-stability  code  LSH,  and  the  linear/nonlinear  PSE  code,  PSH.  Physicad  inter¬ 
faces  axe  set  up  to  amalyze  the  Moo  =  8  boundaxy  layer  over  a  blunt  cone  calculated 
by  using  a  thin-layer  Navier  Stokes  (TLNS)  code,  and  the  flow  over  a  sharp  cone 
at  angle  of  attack  calculated  using  the  AFWAL  Pau’abolized  Navier-Stokes  (PNS) 
code.  While  stability  and  transition  studies  at  high  speeds  are  fair  from  routine,  the 
method  developed  here  is  the  best  tool  available  to  reseairch  the  physical  processes  in 
high-speed  boundary  layers. 
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Chapter  1 
Introduction 


Understamding  and  predicting  the  laminar- turbulent  transition  in  flows  over  aerody¬ 
namic  bodies  is  of  great  importance  to  many  state-of-the-art  technologies  currently 
being  pursued.  Because  skin  friction  and  heat  transfer  in  the  two  states  of  motion 
are  very  different,  it  is  crucial  to  provide  the  practicing  engineer  or  the  researcher 
with  the  best  tools  to  analyze  transitional  flow.  These  tools  enable  one  to  model  and 
control  the  flow  and,  hence,  optimally  design  modern  aierospace  vehicles. 

The  breakdown  of  the  laminar  to  turbulent  flow  is  a  complex  process.  This  pro¬ 
cess  involves  disturbances  in  the  environment,  receptivity  mechanisms  that  internalize 
these  disturbances  into  spatially  varying  basic  flow,  and  streamwise  evolution.  Baring 
bypass  transition,  the  initial  evolution  is  strongly  dependent  on  streamwise  linear  sta¬ 
bility  characteristics.  The  final  breakdown  is  the  result  of  nonlinear  growth,  secondary 
instability,  and  nonlinear  interactions  which  create  increasingly  smaller  scales  of  mo¬ 
tion.  If  the  freestreaun  disturbance  environment  is  known,  the  rem£uning  transition 
process  may  be  predicted. 

Current  tools  for  prediction  may  be  roughly  classified  into  the  following  groups: 

1.  Lump- parameter  methods,  based  on  the  local  boundwy  layer  characteristics 
and  empiriced  correlations. 

2.  The  criterion  (Smith  and  Gamberoni  1956,  Van  Ingen  1956)  based  on  the 
amplitude  ratio  obtained  from  the  linear  stability  theory  and  experimentally 
correlated  N  factors. 

3.  Amplitude  methods  (Liepmann  1943,  Mack  1975a,  b)  based  '''t  only  on  linear 
stability,  but  also  on  the  level  2md  spectral  contribution  of  freestream  distur¬ 
bance  (empiriccilly). 

4.  Methods  based  on  averaged  transport  equations,  with  an  empirical/theoreticaJ 
model  for  the  transition  region. 

5.  Direct  Navier-Stokes  (DNS)  simulations  and  large-eddy  simulations  (LES)  based 
on  the  expensive  numerical  solutions. 

The  first  four  groups  neglect  important  streamwise  vMiations  and  nonlinear  mech¬ 
anisms  of  transition  process.  They  require  sensitive  correlations  and  extensive  effort 
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to  tune  parameters  for  efficient  and  accurate  use  in  complex  design  problems.  The 
last  group  of  numerical  simulations  require  excessive  computer  memory  and  process¬ 
ing  time  to  be  feasible  for  engineering  applications.  The  limitations  for  these  methods 
are  even  more  severe  for  high  supersonic  «ind  hypersonic  flows,  where  experimental 
data  are  scarce  and  the  transition  mechanisms  are  more  complex.  Hence,  alternative 
methods  should  be  developed  for  accurate  and  efficient  engineering  prediction  of  the 
transition  process. 

Alternative  prediction  methods  based  on  the  first  principles  while  taking  all  physics 
into  account  would  also  require  characterization  of  the  freestream  environment,  sur¬ 
face  requirement  such  as  roughness,  and  qualitative  and  quantitative  links  between 
these  disturbances  and  the  initial  conditions  for  the  most  dangerous  modes  of  insta¬ 
bility  mechanisms.  These  data  are  largely  unavailable,  and  receptivity  mech«uiisms 
for  supersonic  flows  are  largely  unknown.  However,  the  existing  baise  of  transition 
data  can  be  used  to  characterize  the  initial  data  at  a  given  environment  (altitude). 
A  proper  integration  method  ol  initial  data  shall  connect  this  initial  data  to  the 
transition  process. 

Therefore,  in  response  to  PRDA  No.  90-07-PMRN,  we  proposed  the  develop¬ 
ment  of  a  transition  prediction  method  for  high  speed  flows  btised  on  the  solution  of 
the  Parabolized  Stability  ^^quations  (PSE).  The  PSE  can  be  solved  efficiently  using 
marching  methods,  while  taking  streamwise  variations  of  the  basic  flow  into  account. 
The  PSE  are,  therefore,  much  more  efficient  than  DNS  methods.  It  also  accounts  for 
the  streamwise  variations  of  the  stability  characteristics,  using  a  streamwise  adapting 
wavenumber  for  fast  variations.  The  shape  of  the  function  accounts  for  the  slow  varia¬ 
tions,  and  only  the  second  streamwise  derivatives  are  dropped.  Moreover,  it  accounts 
for  secondary  instability  and  nonlinear  growth  and  interactions,  enabling  it  to  predict 
the  beginning  of  the  breakdown  process. 

The  concept  was  validated  for  incompressible  flows  (Bertolotti  1991)  at  the  time 
of  the  proposal.  The  main  objectives  set  in  the  proposal  have  been  achieved  and  are 
outlined  as  follow: 

1.  Developing  a  new  set  of  compressible  equations  accounting  for  variable  ther¬ 
mal  properties  of  air  and  body  curvature  by  formulating  the  PSE  in  general 
curvilinear  coordinates. 

2.  Capability  to  generate  amplitude  and  growth  rate  curves  with  ease  for  compar¬ 
ison  with  traditional  methods. 

3.  Developing  a  numerical  approach  for  good  resolution  of  critical  layers  and  wall 
viscous  layers  by  replacing  the  spectral  method  by  a  order  Hermitiain  finite 
difference  method.  Better  distribution  of  grid  points  is  achieved,  and  the  speed 
of  calculations  is  improved  by  factors. 

4.  Setting  the  code  in  a  modular  form  to  permit  use  with  different  types  of  basic 
flows.  This  feature  is  incorporated  by  modular  coding  involving  generic  modules 
which  apply  to  any  stability  problems  and  physics  modules  which  interfaces 
between  different  flows/stability  problems  to  9  general  modules.  In  particulaur, 
physics  interfaces  are  already  in  place  for  flat  plate  flow,  rotating  disk  flow,  flow 
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over  blunt  cones  and  sharp  cones.  The  user  may  add  or  modify  these  modules 
as  new  problems  aue  encountered. 

5.  Addition  of  aJl  nonlinear  terms  to  the  PSE  eind  algorithms  for  introduction, 
dropping  and  interaction  of  different  modes.  Implementation  of  this  module 
enables  the  code  to  predict  secondary  instabilities,  nonlinear  interactions,  and 
beginning  of  brealcdown. 

6.  Ability  to  predict  transition  for  three  dimensioned  (3D)  disturbances  in  two 
dimensional  (2D),  Axisymmetric  or  quasi-3D  beisic  flows.  For  example  linear 
and  nonlinear  evolution  of  disturbances  on  swept  wings  (Stuckert  et  al.  1993) 
eind  disturbance  evolution  2dong  a  suitable  path  on  a  cone  at  angle  of  attack 
(Section  4)  can  be  calculated. 

Two  codes  are  developed  for  this  contract:  a  code  for  local  linear  stability  analysis 
named  LSH  and  a  code  for  analysis  using  the  PSE.  These  two  codes  are  closely 
integrated,  and  use  mostly  the  same  general  modules  and  physics  modules.  Solutions 
from  LSH  can  be  applied  directly  as  initial  conditions  for  PSH  runs.  In  essence,  we 
have  developed  a  pair  of  highly  flexible  codes  for  stability  and  transition  analysis 
that  Ctin  be  adapted  to  a  variety  of  problems  by  the  user.  The  code  can  be  used  with 
ease  for  problems  with  existing  interfaces  (provided  with  the  codes  by  DynaFlow). 
For  new  problems,  one  must  supply  interface  modules  for  metric  terms  of  stability 
coordinate,  basic  state  solution  and  boundary  conditions.  Sometimes  the  user  may  use 
the  existing  metric  file,  but  only  a  new  basic  state  file.  The  tradeoff  for  this  flexibility 
is  that  the  code  requires  the  user  to  have  knowledge  of  stability  and  transition  to  set 
up  a  new  problem.  This  tradeoff  should  not  be  considered  a  disadvantage  because 
this  knowledge  is  helpful  even  for  simple  routine  stability  analyses. 

A  detailed  description  of  the  PSE  formulation  is  given  in  the  next  chapter  followed 
by  the  description  and  results  of  stability  analyses  using  LSH  euad  and  PSH  for  beisic 
flows  for  which  interfaces  are  already  set  up.  Linear  stability  analysis  of  compressible 
flow  over  a  flat  plate  are  presented  in  Chapter  3.  Results  and  discussions  on  linear 
stability  and  PSE  analyses  of  hypersonic  flow  over  a  sharp  cone  at  angle  of  attack 
are  given  in  Chapter  4.  In  Chapter  5,  linear  stability,  linear  PSE  and  nonlinear  PSE 
analyses  of  hypersonic  flow  over  a  blunt  cone  at  zero  angle  of  attack  are  described. 
Some  general  comments  on  the  codes  and  results  axe  given  in  the  last  chapter. 
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Chapter  2 

PSE  Methodology 


The  PSE  predict  the  evolution  of  a  rapidly  changing,  smadl  scale  disturbance  super¬ 
posed  upon  a  steady,  slowly  varying  laminar  flow.  The  undisturbed  lamineir  flow, 
referred  to  as  the  basic  state,  is  presumed  to  be  known,  having  been  computed  ils  an 
exact  or  approximate  solution  of  the  steady-state  Navier-Stokes  equations.  The  total 
flow: 

Q  =  Q  +  Q  (2.1) 

-  the  superposition  of  the  basic  state  and  the  disturbance  -  is,  on  the  other  heind,  to  be 
determined  by  the  solution  of  the  unsteady  Navier-Stokes  equations.  The  unsteady 
Navier-Stokes  equations  really  govern  the  evolution  of  the  disturbance  because  the 
basic  state  is  presumed  known.  Hence,  when  formulated  directly  in  terms  of  the 
disturbances  using  equation  (2.1),  they  are  known  as  the  disturbance  equations. 

The  disturbance  equations  are  nonlinear  because  the  Navier-Stokes  equations  are 
nonlinear.  Although  the  PSE  can  account  for  this  nonlinearity,  it  is  easier  to  first 
consider  the  situation  when  the  disturbance  has  a  vanishingly  small  amplitude.  In  this 
case,  the  equations  are  linear  and  the  disturbeince  coefficients  eu’e  independent  of  time 
because  the  basic  state  is  steady.  Assuming  that  the  disturbance  boundary  conditions 
are  also  separable  in  time,  a  Fourier  transform  can  be  applied  to  the  disturbance  and 
its  equations.  More  generally,  a  discrete  or  continuous  Fourier  tremsform  may  be 
applied  in  each  dimension  in  which  the  basic  flow  and  geometry  do  not  vary.  This 
decomposes  the  disturbance  into  a  multidimensional  spectrum  of  waves: 

0  =  E  (<?m  +  QrA  (2.2) 

TTl— 0 

each  wave  being  governed  by  the  line^u■  disturbance  equations  transformed  from  the 
time  and  space  domain  to  the  frequency  and  wavenumber  domain. 

One  dimension  in  which  a  Fourier  trzmsformation  is  not  strictly  valid  is  the  dis¬ 
tance,  say  s,  in  the  streaimwise  direction.  However,  since  the  basic  flow  veiries  slowly 
with  s,  the  functioned  dependence  of  each  wave  on  s  can  be  decomposed  into  a  number 
of  amplitude  emd  wavenumber-modulated  sinusoids,  the  modulations  varying  slowly 
in  space.  The  governing  equations  for  the  modulated  waves  are  then  obtained  by  per¬ 
forming  yet  another  transformation  of  the  linear  disturbance  equations  to  this  hybrid 
s-space,  s-wavenumber  domain. 
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This  decomposition  is  one  of  the  key  elements  in  the  development  of  the  PSE.  It 
provides  one  with  a  meaningful  way  to  neglect  the  second  derivatives  with  respect  to 
s  of  the  modulated  amplitude  eind  obtain  a  system  of  parabolized  equations  which 
can  be  solved  very  efficiently  using  a  space  marching  technique.  Note  that  usually 
all  of  the  derivatives  of  the  modulations  with  respect  to  s  are  neglected,  yielding  a 
system  of  ordinary  equations  and  an  eigenvalue  problem  to  solve  in  the  transformed 
domain.  However,  for  crossflow  vortices  and  other  longitudinal  vortices,  this  latter 
approximation  can  be  inaccurate. 

The  derivation  of  the  PSE  is  discussed  in  greater  detail  below.  First,  the  pertinent 
variables  are  nondimensionalized.  The  Navier-Stokes  equations  are  then  presented  in 
an  inertial  Cartesian  coordinate  system  and  linearized  around  the  laminar  basic  state 
to  yield  the  linear  disturbance  equations.  A  coordinate  transformation  is  then  made 
to  identify  meaningful  independent  variables  for  the  Fourier  tremsformations  and  sub¬ 
sequent  approximations.  The  coordinate  transformation  also  enables  the  boundary 
conditions  to  be  enforced  easily  and  accurately.  Finally,  some  of  the  details  of  the 
nonlinear  analyses  are  also  presented. 


2.1  Nondimensionalizations 


X 

= 

(2.3) 

t 

= 

(2.4) 

= 

(2.5) 

p 

= 

Plip-ooU^J) 

(2.6) 

T 

= 

rir^ 

(2.7) 

P 

= 

pVpI. 

(2.8) 

h 

= 

hyp-T:, 

(2.9) 

P 

= 

pVpI. 

(2.10) 

X 

= 

y/pl. 

(2.11) 

K 

= 

(2.12) 

= 

(2.13) 

= 

(2.14) 

Note  that  the  Prandtl  number  Pr  =  ^*c*/k*  =  fiCpIn. 

For  the  cone  calculations  considered  in  this  report,  the  reference  temperature  is  the 
free-stream  static  temperature;  the  reference  viscosity  is  the  free-stream  viscosity;  the 
reference  density  is  the  free-stream  density;  the  reference  velocity  is  the  free-stream 
velocity:  and  the  reference  length  is  the  length  of  the  model  divided  by  the  square 
root  of  the  Reynolds  number  based  upon  the  reference  density,  velocity,  viscosity,  and 
model  length. 

With  these  nondimensionalizations,  the  Navier-Stokes  equations  in  an  inerti£il 
Cartesian  coordinate  system  can  be  written  as  follows. 
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2.2 


Navier- Stokes  Equations 


dp  dpu‘‘ 

du^  dp  1  95^* 

^  dt  ^  dx^  dx^  Re  dx^ 

(2.15) 

(2.16) 

9k  US  (dp  e9?\ 

RrT^\dt'^  ai‘j 

_  -ia,‘  1  us 

Redx'‘  ReRrT;^ 

(2.17) 

(For  an  ideal  gas,  the  speed  of  sound  a*  =  and  U^fR!‘T^  =  Here, 

the  index  k  6  (1,2,3)  and  x=  {x^,x^,x^)  is  the  position  vector  in  the  inertial  Carte¬ 
sian  coordinate  system.  The  viscous  stress  tensor  5^*  is 

\dx^  dx^J  dx^ 

(2.18) 

The  components  of  the  heat  flux  axe 

it  ^ 

(2.19) 

Finally,  the  viscous  dissipation  is 

$  -5>‘— 

’  ^  ax'- 

(2.20) 

2.3  Disturbance  Equations 

The  disturbance  equations  are  obt2dned  by  substituting  the  disturbed  flow  field  (2.1) 
into  the  Navier-Stokes  equations.  Nonlinear  terms  are  then  expanded  in  a  truncated 
Taylor  series  about  the  basic  state  so  that  each  is  in  the  form  of  a  sum  of  products 
of  the  disturbance  variables. 

A  typical  term  in  the  momentum  equations  is  of  the  form: 
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The  nonlineax  disturbance  equations  are  obtained  by  similarly  expanding  all  of  the 
terms  in  the  Navier-Stokes  equations.  The  equations  will  then  contain  three  different 
types  of  products.  Some  products  will  involve  basic  state  quantities  only  (i.e.,  the 
first  of  the  terms  in  equation  (2.22)),  some  will  be  linear  in  the  disturbances  (i.e., 
the  second  and  third  terms  in  equation  (2.22)),  and  the  rest  will  be  nonlinear  in  the 
disturbances  (i.e.,  the  last  term  in  equation  (2.22)).  These  terms  in  the  equations  can 
then  be  grouped  accordingly. 

In  each  governing  equation  the  group  of  basic  state  terms  will  sum  to  zero  because 
the  basic  state  is  tissumed  to  be  an  exact  solution  of  the  Navier-Stokes  equations. 
The  linear  terms  can  then  be  placed  on  the  right  hand  sides  of  the  equations,  and 
the  nonlinear  terms  on  the  left  hand  sides  of  the  equations.  The  nonlinear  terms 
are  then  expanded  in  Fourier  series  using  the  elementary  modes  as  the  basis 
functions.  Orthogonality  conditions  between  the  basis  functions  (or  the  principle  of 
harmonic  balancing  -  see  Herbert  (1991)  and  Bertolotti  (1991)  for  more  details)  are 
then  enforced  to  obtain  the  equations  governing  the  evolution  of  each  mode. 

In  the  nonlinear  case,  the  equations  for  different  modes  may  be  coupled  to  each 
other.  However,  the  linear  terms  in  the  equations  are  the  same  regardless  of  the 
mode  amplitudes.  Hence,  the  linear  terms  are  presented  below  in  the  form  of  the 
linear  disturbance  equations. 

Denote  the  linearized  material  derivative  acting  on  any  quantity  by: 

DQi  _aQ‘  u»Q‘  ,dQi 

Dt  ~  dt  ai* 

where  Q  =  =  (p,T',u^^^^u^). 

The  linear  disturbance  equations  in  inertial  Cartesian  coordinates  are: 


Dp  du’’ 


(2.25) 


d^h  DQ^\.,  dh  DQ^' 

Dt  Dt 


,  dh  DQ^ 

+  pI. 


U-J_  g  dp  DQ^ 


f^^dQ^  Dt  BrT^^^dQ^  Dt 

^  ^  d^p 

s.T’*  2-^  \  2^ 


+ 


,Du^ 

^~Di 


^  Dt 


Vrti 

1  dq'^ 

Re  dx’’ 


+ 


Du^  .  ^  d^p  dQ^\  .. 

^  Dt  dQ^dQ^  dxi )  ^ 


^  ^  dp  dQ^ 


1  dS^’^ 


^  dQ"^  dxi  Re  dx^ 


=  0 


(2.26) 


(2.27) 
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The  linear  disturbance  in  the  divergence  of  the  heat  flux  is: 


dQ^dT\., 
ax*  dQ‘dQ^  ax*  ax* )  ^ 

^  dK  dQ'  dT  a/C  dQ^  df 

^  ^  a(p’’  ax*  ax*  ^  ^  a^*^  ax*  ax* 

^  at?’"  /  ax*ax*  ^  ax*ax* 

The  linear  disturbance  in  the  divergence  of  the  shear  stress  tensor  is: 


(2.28) 


a5J* 

ax* 


+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 


yfy  ay  ag- 

^  Utt  ag*ac'  ax*  vax*  ax^  y  /  ^ 

iiXk'i  dQ‘dQ^  ax*  J  ^ 

^  dQ^  ax*  \ax*  ^  axJ  j 

a^’-  ax*  vaxJ ) 

^  a//  aQ^  /au^  ati*\ 

dQ'  ax*  vax*  ax>  j 

^  a^’' ax*  Vaxj  j 


a/i  a^u*  a^u* 

aQ*"  dx^dxj  ^  dx^dx^ 


dx^dx^ 


a^u* 

dx^dx^ 


(2.29) 


Next,  the  coordinate  system  is  transformed  to  enable  the  accurate  implementation 
of  the  PSE  approximation  and  boundary  conditions. 


2.4  Transformation  of  the  Independent  and  De¬ 
pendent  Variables 


Let  the  coordinate  trsmsformation  be  given  by: 

=  Am 

t  =  T 


(2.30) 

(2.31) 
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where  the  body  surface  is  =  0.  The  index  k  €  (1,2,3)  and  *=  i^)  is  the 

position  vector  in  the  inertial  Cartesian  coordinate  system. 

The  derivatives  of  the  disturbance  variables  with  respect  to  can  easily  be  com¬ 
puted  using  the  chain  rule  for  partial  differentiation.  Before  this  is  actually  done, 
however,  the  dependent  vau’iables  are  also  transformed.  This  is  because  the  approxi¬ 
mations  used  to  solve  for  the  disturbances  are  not  valid  for  every  possible  choice  of 
dependent  variables  that  one  may  choose  to  work  with.  In  the  actual  software  that 
we  have  developed,  a  number  of  reasonable  options  have  been  implemented.  One  that 
is  commonly  chosen  uses  the  density  and  temperature: 

Qi  =  f^itr)exp{ier^ie,e.r)}  (fc  =  l,2)  (2.32) 

and  the  physical  values  of  the  contravariant  velocity  components  in  the  curvilinear 
coordinate  system  (i.e.,  Borisenko  and  Tarapov  [1979],  pp.  23-35): 

1 

Qi  =  (2.33) 

{k  =  3,4,5) 


Here,  is  the  jth  component  of  the  shape  function  f  for  the  elementary  mode. 
The  repeated  index  j  in  the  previous  equation  indicates  a  summation  from  1  to  3.  gjj 
(no  summation)  is  the  diagonal  element  in  row  j  of  the  covariant  metric  tensor  gif 


~'dfW 


(2.34) 


where  the  repeated  index  k  also  indicates  a  summation  from  1  to  3.  0m  is  the  (com¬ 
plex)  phaise  for  the  mode  and  is  assumed  to  be  a  function  of  ^^(*,t), 
and  r,  but  not  of  ^^(®,t).  Its  gradient  is: 

d9 


(2.37) 

This  assumed  form  for  the  disturbance  modes  is  fundamental  to  the  Parabolized 
Stability  Equations  and  has  arisen  as  an  extension  of  the  local  linear  stability  analyses. 
The  local  linear  stability  analyses  are  more  restrictive  because  they  formally  require 
that  the  disturbance  equations  be  separable  in  t,  and  Since  the  PSE  account 
for  the  variations  in  the  shape  functions  of  the  different  modes,  they  are  not  bound 
by  this  restriction.  Moreover,  they  can  also  account  for  inhomogeneous  boundau’y 
conditions  and  nonlinearity. 


de 

_d9 

dr 


2.5  Nonlinear  Terms  and  Modal  Decomposition 

The  means  by  which  the  nonlinear  terms  in  the  disturbance  equations  are  incorporated 
can  now  be  explained  more  fully.  Assume  for  simplicity  that  Cartesian  coordinates 
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are  used.  Substituting  equations  (2.32)  and  (2.33)  into  equation  (2.23),  one  obtains 
sums  of  products  of  the  form: 


[Prn  +  Pm)  H  ^ 

msO  /=0  t  ^ 


nmnm  ( _  g-j  _  g~J  ^  g.j  ^  g^J 


=  EE 

m=0/=0  \  ^ 


dt  dt 


(2.38) 


If  the  underlying  laminar  flow  is  stationary  and  invariant  in  the  circumferential  direc¬ 
tion,  then  u>m  aJid  /3m  are  constants  for  e2u:h  m.  Moreover,  if  one  chooses  an  equally 
spaced  “grid”  in  the  frequency  and  circumferential  wavenumber  domains,  then  the 
frequency  and  circumferential  wavenumber  for  each  mode  are  integral  multiples  of 
some  corresponding  fundamentals  Aw  and  A/3,  respectively: 


=  j(m)Aw  (2.39) 

/3m  =  k(m)A/3  (2.40) 

The  product  of  two  typical  modes  mi  and  m2  in  equation  (2.38)  above  then  yields 
terms  proportional  to 


exp  [i  {;(mi)  -f  i(m2)}  (Aw)t  -1-  i  {A:(mi)  -t-  ^(mj)}  {A0)f]  (2.41) 

Due  to  the  orthogonality  of  the  sinusoidal  functions,  this  creates  a  nonlinear  forcing 
term  for  the  mode  m3  which  simultaneously  satisfies 

i(m3)  =  j(mi)-b;(m2)  (2.42) 

k{m3)  =  k{mi)  fc(m2)  (2.43) 


Products  of  one  complex  mode  with  the  conjugate  of  another  mode  generate  differ¬ 
ences  in  the  phase  functions  instead  of  sums,  etc.  Also,  higher  order  (i.e.,  cubic) 
nonline2irities  axe  treated  in  the  same  way,  yielding  the  nonlinecir  equations  governing 
the  evolution  of  each  mode. 

The  final  governing  equations  in  the  curvilinear  coordinate  system  axe  obtained 
by  substituting  equations  (2.32)  and  (2.33)  into  the  disturbance  equations  and  trans¬ 
forming  the  derivatives  using  the  chain  rule  for  paxtiid  differentiation.  The  nonlinear 
terms  are  computed  as  described  above,  and  then  the  PSE  approximation  is  made. 

2.6  PSE  Approximation 

In  linear  stability  theory  the  variations  of  the  shape  function  fm  and  otm,  0m,  and  Um 
with  and  axe  usually  neglected  in  order  to  reduce  the  governing  system  of  peutial 
differential  equations  to  a  system  of  ordinary  differential  equations.  The  boimdaxy 
conditions  are  also  usually  assumed  to  be  homogeneous.  Two  of  the  three  parameters 
am,  0m  and  Um  must  then  be  specified  and  a  complex  eigenvalue  problem  solved  for 
the  third. 
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The  PSE  methodology  still  requires  the  specification  of  constant  u>m  and 
However,  the  variations  of  the  shape  function  with  respect  to  are  not  neglected. 
Hence,  the  partition  of  the  vaxiations  of  the  disturbance  with  respect  to  between 
the  shape  function  and  the  exponential  factor  exp{i^TO(^S^^,'r)}  must  be  specified. 
A  method  for  doing  this  is  described  by  Herbert  (1991).  It  is  motivated  by  a  desire 
to  make  “small”  relative  to  other  terms  in  the  governing  equations  so 

that  it  may  be  neglected. 

Consider  any  component  /  of  the  m**  mode.  Its  logarithmic  derivative  is: 

^  {M/m(^)exp{^0,„(^^^^T)}]} 

=  +  (2.44) 

Since  depends  upon  both  and  this  equation  in  general  implies  a  dependence 
of  a  upon  both  I  and  4^.  If  one  specifies  a  particular  value  of  1,  the  dependence 
of  a  upon  may  be  eliminated  by  multipl3ang  equation  (2.44)  by  a  weight  factor 
I/4P  =  fLPm  and  integrating  over  the  entire  domain  of 


Sf""  l/il^afr  {ln[/^(^)exp(z^,n)]}  df 


/o'”"  ifLi^de 

(2.45) 

The  partition  is  now  specified  by  setting 

filas  -t  dfL  , 

L 

(2.46) 

Equation  (2.45)  then  becomes: 

l/inP^r  {lnl/i(«)exp(i9„)l}<ie^  , 

- - - -  =  lOjn 

/o"" 

(2.47) 

and  is  used  to  specify  a(^^). 

Note  that  the  partition  depends  upon  /  aind  the  grid  used  for  the  calculations. 
However,  different  choices  of  I  and  different  grids  should  lead  to  the  same  physical 
solution  Q  to  within  the  order  of  error  incurred  by  the  PSE  approximation.^  This 
approximation  is  that  the  derivative  is  negligibly  small  in  comparison 

with  the  remaining  terms  in  the  disturbance  equations  and,  hence,  can  be  dropped.^ 
The  linear  PSE,  thus,  can  be  obtained  from  the  disturbance  equations  by  neglecting 
the  (viscous)  terms  involving 

For  any  particular  disturbance  calculation  the  equation  of  state,  transport  prop¬ 
erties,  boundary  conditions  and  initial  conditions,  grid  =  ^*(®,  <),  and  basic  state 
must  be  specified.  These  are  described  next. 

^We  have  used  /  =  3. 

^Note  that  the  shape  function  is  independent  of  for  aui  axisymmetric  body  at  zero  angle  of 
attack  because  both  the  geometry  and  flow  field  are. 
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2.7  Constitutive  Equations 

The  gas  is  cissumed  to  be  air  ajid  is  modeled  as  perfect  in  this  analysis.  Hence,  the 
equation  of  state  is: 


P{P,T)  = 


The  enthalpy  h  is; 


h{T)  = 

where  7  =  1.4.  Sutherland’s  law  is  used  for  the  dynamic  viscosity: 

73/2 


(2.48) 

iT 

7  -  1 

(2.49) 

7 

7-1 

(2.50) 

/i(T)  =  (l  +  nOAK/T^ 


J  +  mAK/T^ 

The  second  coefficient  of  viscosity  is  obtained  from  Stokes’  hypothesis: 


A  =  -3. 


(2.51) 


(2.52) 


Finally,  the  thermal  conductivity  is  computed  assuming  that  the  Prandtl  number  is 
a  constant  equal  to  0.72: 


k  =  ^ 
0.72 


(2.53) 


Other  constitutive  equations  for  the  thermophysicid  properties  of  air  cam  adso  be 
selected.  See  the  Appendix:  Constitutive  Equations. 


2.8  Boundary  Conditions  and  Initial  Conditions 

The  following  general  conditions  are  applied  in  boundary .  ins  files  provided  with  the 
built-in  sample  applications.  At  the  wall,  disturbance  temperature  and  velocities 
are  set  to  zero.  At  farfield,  disturbance  temperature  Q^,  velocities  in  and 

directions  are  set  to  zero.  Then  ^^-gradient  of  disturbance  velocity  in  direction 
is  determined  from  the  continuity  equation  aissuming  that  density  and  gradients 
of  density  are  zero  implicitly.^  No  explicit  boundary  condition  on  density  is  required 
at  either  end  of  the  boundauries.  For  the  current  version  of  the  code,  asymptotic 
conditions  axe  implemented  only  for  blunt-cone  calculations  using  orthogonal  body- 
intrinsic  coordinates.  Boundary  conditions  for  the  stationary  {u  =  0)  modes  may 
be  different  from  those  of  oscillatory  disturbances.  For  example,  when  the  basic 

*The  user  may  choose  to  set  this  velocity  component  to  zero  for  oscillatory  modes  with  little  or 
no  change  in  the  solution. 
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l2iminax  flow  satisfies  the  adiabatic  wall  boundary  condition,  normal  gradient  of  T  for 
stationaiy  modes  may  be  set  to  zero  (see  Mack  1984). 

The  disturbance  velocity  and  temperature  gradient  are  <issumed  to  vanish  at  the 
surface.  Initial  conditions  are  obtained  by  making  the  locally  queisi-paxallel  flow 
approximation:  the  shape  function  is  assumed  to  depend  only  upon  the  covariant 
components  and  0rn  of  tbe  wavenumber  vector  are  adl  assumed  to  be  constant; 
and  all  of  the  coefficients  in  the  disturbance  equations  which  depend  on  and  axe 
evaluated  at  the  location  of  the  initial  data. 


2.9  Numerical  Method 

The  PSE  approximation  enables  a  solution  for  the  disturbance  to  be  obtained  effi¬ 
ciently  by  integrating  the  equations  in  the  direction  with  a  space  m«irching  cilgo- 
rithm.  Here,  the  implicit  backward  Euler  method  is  used  for  discretization  of  deriva¬ 
tives  with  respect  to  The  nonlinear  terms,  however,  have  been  treated  explicitly  as 
this  saves  a  significant  amount  of  computational  time.  Also,  our  earlier  calculations 
had  used  a  spectral  collocation  method  for  discretization  of  derivatives  with  respect 
to  We  have  since  begun  to  use  a  seventh  order  accurate  compact  finite  difference 
scheme  on  the  interior  which  leads  to  a  system  of  block  tri diagonal  equations.  This 
system  of  block  tridiagonal  equations  can  be  solved  much  more  efficiently  than  the 
full  matrices  obtained  with  the  spectral  method  while  the  seventh  order  accuracy 
maintains  the  quality  of  the  solution.  In  summary,  the  solution  algorithm  is  now 
much  more  efficient. 

Nevertheless,  since  even  the  equations  for  a  single  mode  are  nonlinear  in  due 
to  its  appeareince  in  the  differential  operators  acting  upon  the  shape  function,  Om 
must  first  be  guessed  at  each  new  streamwise  station,  updated  by  imposing  the  norm 
on  /,  and  then  iterated  upon  until  convergence.  These  iterations  typically  converge 
rapidly,  but  the  rate  of  convergence  depends  upon  the  step  size  and  flow  and,  in  the 
case  of  nonlinear  analyses,  on  the  amplitudes  of  the  disturbances. 
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Chapter  3 


Compressible  Bound2iry  Layer 
over  the  Flat  Plate 

The  stability  of  the  compressible  boundary  layer  on  a  flat  plate  is  investigated  as  a 
test  case  for  the  insert  files  which  incorporate  general  curvilinear  coordinates.  Two 
different  coordinate  systems  aje  used.  The  fiist  of  these  is  a  stretched  rectilinear  grid 
whose  point  distribution  matches  that  of  the  wall  normad  similarity  variable  adong  the 
inflow  boimdary.  The  second  coordinate  system  is  curvilinear  and  nonorthogonal.  It 
uses  the  wadi  normad  similarity  variable  and  the  distance  measured  from  the  leaiding 
edge  of  the  plate  as  its  two  coordinates. 

While  we  now  use  the  Hermitian  finite  difference  method  to  discretize  and  solve 
the  disturbance  equations,  the  results  presented  in  this  section  were  obtained  with 
the  slower  but  extremely  accurate  spectrad  collocation  technique.  The  high  accuracy 
of  the  spectrad  method  enables  very  careful  comparisons  to  be  made  which  can  reveal 
errors  in  the  formulation  or  software. 

In  addition,  the  linear  disturbance  equations  are  formulated  in  terms  of  (p,  T, 
u,  V,  xv),  where  p  is  the  disturbance  pressure,  T  the  disturbance  temperature,  and 
(u,  u,  u;)  the  disturbance  velocity  vector  in  Cartesian  coordinates.  This  selection  of 
dependent  variables  is  useful  because  it  makes  it  easier  to  study  incompressible  flows 
simply  by  changing  the  equation  of  state  from  the  ideal  gas  law  to  p  =  1.^  However, 
the  subsequent  nonlinear  work  uses  as  dependent  variables:  (p,T',ti,t;,u;),  the  same 
set  as  appears  in  the  formulation  of  the  disturbance  equations  in  Chapter  2.  This 
latter  choice  makes  it  easier  to  formulate  the  nonlinear  terms  because  it  gives  rise 
to  products  of  sums  of  the  basic  state  and  disturbance  modes  rather  than  quotients. 
The  drawback  is  that  it  seems  to  require  many  more  Chebyshev  polynomials  than 
does  the  pressure  formulation  for  the  same  level  of  accuracy.  We  observed  this  when 
we  confirmed  the  changes  in  the  linear  equations  by  comparing  the  results  from  the 
density  formulation  to  those  from  the  pressure  formulation. 

*  Incompressible  flows  were  studied  and  tested  prior  to  performing  the  compressible  flow  tests 
outlined  here. 
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3.1  Basic  State  Similarity  Solution 

The  basic  state  compressible  boundary  layer  is  specified  to  be  the  Blaisius  similarity 
solution.  This  section  presents  some  of  the  details  of  the  similtirity  transformation 
£ind  resulting  equations. 


3.1.1  Governing  Equations  and  Similarity  Transformation 


The  two-dimensional  compressible  boundary  layer  equations  for  steady  flow  over  a 
flat  plate  are: 


dpu  dpv 
dx  dy 


du 


du 


pu—  +  pv-^ 
ox  ay 

dp 

dy 


dh  dh 


=  0 

?}i] 

Reidy  \  dy) 

=  0 

1  d  (.  dT\ 


(3.1) 

(3.2) 

(3.3) 

(3.4) 


(See,  for  instance,  White  [1974].) 

Here  x  is  the  distance  from  the  leading  edge  of  the  plate  and  y  ij  the  distance 
measured  normal  to  it.  u  and  v  are  the  corresponding  velocity  components,  p  is  the 
density  of  the  g«is,  and  h  is  its  enthalpy,  p  is  the  static  pressure  and  T  is  the  static 
temperature,  p  cind  k  are  the  dynamic  viscosity  and  thermal  conductivity  of  the  gas, 
respectively. 

The  gas  is  assumed  to  be  ideed  so  that  the  enthalpy  is  a  function  only  of  the 
temperature.  The  transport  properties  are  eilso  assumed  to  be  functions  of  the  tem¬ 
perature  only.  _ 

All  length  scales  have  been  nondimensiondized  by  =  yj XqPI! PlU*  where  is 
the  nominal  distance  from  the  leading  edge  of  the  plate  and  pi,  pl,  and  f7*  are  the  dy¬ 
namic  viscosity,  g2is  density,  and  flow  speed,  respectively,  at  the  boundary  layer  edge. 
The  pressure  has  been  nondimensionalized  by  plU*^,  but  the  temperature  has  been 
nondimensionalized  by  its  edge  V2due  T*  and  the  enthzdpy  hzis  been  nondimensional¬ 
ized  by  R’Tg.  (/?*  is  the  gas  constant.)  The  viscosity  heis  been  nondimensionalized 
by  pI  and  the  thermal  conductivity  by  plR".  (Note,  then,  that  ihe  Prandtl  num¬ 
ber  =  p^c^jk*  =  pCpfk,  where  Cp  =  dh/dT  is  the  specific  heat  at  constant  pressure 
nondimensionalized  by  the  gais  constant  if.) 

Finally,  Rti  =  plU’LUpl  is  the  Reynolds  number,  7e  =  Cpel{cj,e  —  1)  is  the  ratio 
of  specific  heats  at  the  boundary  layer  edge,  and  Me  =  is  the  boundary 

layer  edge  Mach  number. 

A  similarity  solution  for  these  equations  can  be  found  (again,  see  White  [1974]) 
subject  to  the  boundary  conditions: 


u(0,y)  =  1 


(3.5) 
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r(o,y)  =  1 

u(x,0)  =  0 
v(i,0)  =  0 
dT 

^(T.O)  =  0 

limu(x,«)  =  1 

y-^OO  '  ^  ' 

lim  T{x,y)  =  1 

y-*00  *'  ' 


(3.6) 

(3.7) 

(3.8) 

(3.9) 

(3.10) 

(3.11) 


This  involves  (a  modification  of)  the  Illingsworth  trainsformation  of  the  independent 
variables: 


X  =  ^ 

y  =  I  T{s)ds 
V  xq  Jo 


(3.12) 

(3.13) 


Here,  T  =  T{t))  is  assumed  to  be  a  function  only  of  rj  cind  not  of  x.  Also,  make  note 
of  the  fact  that  tq  =  Rcl  with  the  present  choice  of  length  scale. 

The  derivative  matrix  of  x  and  y  with  respect  to  ^  amd  r)  is 


and  its  inverse  is: 


dx 

dx 

a? 

dn 

ay 

a{ 

dv 

§L 

dx 

dy 

§IL 

in 

dx 

dy 

1  0 


(3.14) 


^  n 


f^T{s)ds 


/2xfxoT(rt) 


(3.15) 


Now  introduce  the  stream  fimction  rp: 


In  terms  of  ^  and  77, 


1 

\I2xJxqT 
di  ^  2xT . 


d(  2xT  Jo  dv 


(3.16) 

(3.17) 

(3.18) 

(3.19) 


Since  the  pressure  is  constant  through  the  boundary  layer  amd  the  g«is  is  ideal, 
p{x,y)T{x,y)  —  1.  Hence,  equation  (3.18)  can  be  written: 


1  dip 
\j2xlxo 


(3.20) 
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Assuming  that  u  is  a  function  df/di;  of  -q  alone  (where  dfjdq  =  0  at  77  =  0  to  satisfy 
the  no  slip  condition), 


yj2i  jxo 

Integrating  with  respect  to  q  gives 

Substituting  this  into  equation  (3.19)  and  using  pT  ~  1  gives 


(3.21) 


(3.22) 


“  "  s/f  lo  "  fMTM)  -  mTM  (3.23) 

Since  there  is  no  wall  suction,  one  must  have  =  0  in  order  to  satisfy  i;  =  0  at 
the  wall.  For  simplicity,  then,  set  ^(^)  =  0. 

To  summarize,  the  stream  function  is: 


Xq 


The  tangential  velocity  component  is: 


“ = !<''> 

and  the  normal  velocity  component  is: 

Using  these  expressions,  the  material  derivative  operator  can  be  written 

D  d  d 

Dt  ''dx'^^dy 

-  <’'>3? '  ITa^ 

The  boundary  layer  equations  thus  become  (with  p  —  IjT) 


-IL  =  /  1  U  1  \d  (pf" 

2xT  VRei)  \2Txlxol  dq 
fc,T 


(¥) 


Sf  ■  (5:)(iik){i(¥)*'*<(f)in’} 

Simplifying  (while  using  lo  =  Rei), 


(3.24) 

(3.25) 

(3.26) 

(3.27) 

(3.28) 

(3.29) 

(3.30) 

(3.31) 

(3.32) 
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If  the  specific  heat  and  the  Prandtl  number  are  constant,  then  Cp  =  7e/(7e  —  1)  and 
k  =  /iCp/Pr  so  that 

(f-) (f) 

These  equations  were  solved  using  the  Chebyshev  spectral  collocation  technique 
in  the  bstate .  ins  file  for  the  stability  analysis. 


3.2  Stability  Analysis 

This  section  describes  some  stability  analyses  of  the  flat  plate  compressible  boundary 
layer.  Two  different  coordinate  systems  are  used  in  conjunction  with  three  different 
sets  of  dependent  variables.  Parallel,  quasi-parallel,  and  PSE  ciilculations  aie  done 
for  most  of  the  cases.  Finally,  some  nonlinear  PSE  computations  eu-e  compared  with 
the  results  of  a  linear  secondary  stability  analysis  (El  Hady  [1991]). 


3.2.1  Coordinate  Systems 

Both  a  rectilinear  grid  and  one  obtjuned  from  the  similarity  transformation  are  used. 
Since  the  rectilinear  grid  makes  use  of  some  of  the  similarity  transformation  relations, 
the  similarity  tr£msformation  is  presented  first. 


Similarity  Transformation 

The  similarity  transformation  has  already  been  described  for  the  most  part  (equations 
(3.13),  (3.14),  and  (3.15)).  These  equations  are  presented  again  here  for  convenience: 


I  =  ^  (3.35) 

y  Xq  Jo 

Again,  T  =  T{ti)  is  assumed  to  be  a  function  only  of  rj  and  not  of  x.  Remember  too 
that  xo  =  Rcl  with  the  present  choice  of  length  scale:  Lq  =  y^Xo/r*/pji7*. 

The  derivative  matrix  of  x  md  y  with  respect  to  ^  and  rj  is 


and  its  inverse  is: 


dx 

dx 

9ti 

djl 

3{ 

dv 

§l 

di 

dx 

dy 

tL 

dx 

dy 

2xrJ0  V  ;  y/2xixonn) 


(3.37) 


(3.38) 
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The  second  derivatives  of  x  =  and  y  —  with  respect  to  (  =  ^^,ri  =  are: 


d^x 

-{ 

d^y 

_ 

- 

•  d^y  1  ^ 

[  i\/fr(,)  7|r(,)  j 

The  third  derivatives  of  y  with  respect  to  ^  are: 

a  [  yy  1  ^  ( ^^/¥,!sn‘)ds  ^^T{v) 

d(  [dm‘\  ^  ^(y)  h^Tiv) 

±  [_^!l1  =  (  ivf  ^'(y)  \ 

dr,  [af ac'J  1  i^r(y)  yf^"(y)  j 


(3.39) 

(3.40) 


(3.41) 

(3.42) 


The  covariant  and  contravariant  components  of  the  metric  tensor  as  well  as  the 
Christoffel  symbols  of  the  first  and  second  kinds  are  given  here  too.  The  deriva¬ 
tives  of  many  of  these  quantities  are  also  given.  They  are  noi  needed  as  input  in 
any  way;  they  are  included  for  compzurison  to  the  values  which  axe  computed  by  the 
program.  In  fact,  the  computed  results  were  compared  with  those  given  here  using 
the  UNIX  debugger  -  they  are  in  complete  agreement. 

The  covariant  components  of  the  metric  tensor  are 


w-'l- [a?  a{>  J  \  ^ffT(s)ds  2(ft)rn-7)  j 

The  contravariant  components  of  the  metric  tensor  are 

f  ‘^1  =  ffl-r'  =  f  ^  2’(5)ds 

V  2xt}v)  T{s)ds  2(x/xo)T^(v)  {.2x^1,)  fo  T{s)ds^ 

The  Jacobian  of  the  coordinate  transformation  is 


(3.43) 


(3.44) 


J  =  y/g  =  vdet[5ij]  = 


Hence, 


0 


(3.45) 

-iln(xo)  +  ll>(r(,)) 

(3.4«) 

IM  \ 

’  T{v)  ) 

(3.47) 

0  \ 

(n,)-*^)  j 

(3.48) 
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The  ChristofFel  symbols  of  the  first  kind  are  defined  by 


[st,p]  = 

dgps  dg,t\ 

2  V  ) 

(3.49) 

They  an'e  symmetric  in  s  and  i.  The  nonzero  elements  for  the  similarity  transformation 
are: 

[11.11 

=  4x4  [L 

(3.50) 

(12. 1] 

(3.51) 

(22, 1| 

=  rns)i. 

xq  Jo 

(3.52) 

(11.21 

=  -J^'>U^Tis)ds 

2xxo  Jo 

(3.53) 

(12.2] 

_  nv) 

Xo 

(3.54) 

(22,2) 

=  2(^)r(,)r(,) 

(3.55) 

The  Christoffel  symbols  of  the  second  kind  are  defined  by 

III 

•3 

Co 

(3.56) 

They  are  also  symmetric  in  s  and  t.  The  nonzero  elements  for  the  similarity  trans¬ 
formation  are: 

f  2  1 
111) 

(3.57) 

f  2  1 

tl2J 

1 

^  “  2x 

(3.58) 

f  2  1 

l22J 

II 

_ .  ^ _ 

(3.59) 

The  derivatives  of  the  ChristofFel  symbols  of  the  second  kind  with  respect  to  ^  and  t/ 
are: 

^  I  1  1  I  "  (  2x3T{v)  ^0  T{s)ds  ,  ^T{v)  “  7^  Jo  T{s)ds^  )(3.60) 

A}  =  ’  “)  <*■'**> 
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Again,  aJl  of  the  above  vadues  are  computed  correctly  by  the  files  used  with  the 
programs  LSH  and  PSX. 


Rectilinear  Coordinate  System 

The  rectilinear  coordinate  system  is  defined  by  the  transformation 

(  =  ^  =  x 

V  =  Tl{xo,y) 

Hence, 


dx  dy 
drj  df) 
dx  dy 


0  |^(®o,y) 


(3.63) 

(3.64) 


(3.65) 


Higher  derivatives  are  computed  similarly  and  axe  not  given  here.  However,  since 
the  basic  state  is  computed  in  terms  of  ^  and  it  is  convenient  to  figure  out  the 
transformation  laws  between  the  4,  V  V  coordinates.  These  transformation 

laws  are  given  below. 

^  =  1  (3.66) 


r  =  0 

dr} 


dr}  dx  drj  dy 

dx  d^  dy  d^ 

drj  dx  dfj  dy 

dx  drj  ^  dy  drj 


or,  using  (3.65), 

/  §l  Si\ 

a(  I  _ 

[  M  §R  J 

\  d(  dr,  / 

The  inverse  of  this  matrix  is: 

/  ^  M  \  / 

di  1  _  * 


|;(xo,!i)t(x,y)  |;{io,y)|;(i,y) 


lit  Hz'  I 

or,  using  (3.37), 


^  14  \ 

9€  dv  \  _ 

^  I?  I  “ 

94  dv  } 


_=1_  jr/  \  j  i  mxo,y)) _ 


TMxo,y 


(3.67) 

(3.68) 

(3.69) 

(3.70) 

(3.71) 

(3.72) 


(3.73) 


2xT(7j)  /o  T{s)ds 


/x/xoT{r,) 


(3.74) 
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The  chain  rule  transformation  law: 

d  d  d  drj 

(3.75) 

d  d  d  drf 

di)  dy  ^  drj  dy 

(3.76) 

is  thus: 

dl  di  2xT(Ti)Jo 

(3.77) 

d  T{r})  d 

drj  ^xIxoT{t))9v 

(3.78) 

The  second  derivative  operators  can  be  shown  to  be 

52 

de 

ST(s)ds(  T'M  ^ 

de  {2xTMf  V  ^  TM  1  *  '  7  a. 

(3.79) 

didrj 

TM)  [S’  \  [T'M  ,1  s 

y/x/xoT{Ti)Wdv  2x[T^{r})Jo  ] 

(3.80) 

2xT{v)Io 

(3.81) 

drj^ 

1  \riT})  1  /r(^)yr(»/)l  d 

^jxlxt,  r{ri)  ^x/xo\S'{^)}  S{ri) ^  dri 

(3.82) 

\  m  Y  S’ 

l^x/xoTM}  Sv’ 

(3.83) 

This  coordinate  system  was  only  used  in  actual  computations  for  the  local  analyses 
(parallel,  quasi-parallel).  Similarity  and  Cartesian  coordinates  are  used  for  the  PSE 
calculations. 


3.2.2  Dependent  Disturbance  State  Variables 

As  mentioned  at  the  beginning  of  this  section,  three  different  sets  of  dependent  vari¬ 
ables  can  be  used. 

The  first  set  is:  (p,  T,  u,  v,  w),  where  p  is  the  pressure,  T  the  temperature,  and  u,  u, 
w  the  velocity  components  in  the  i,  y,  and  z  (spanwise)  coordinate  directions,  respec¬ 
tively.  The  ■  overstrike  denotes  the  disturbance.  These  variables  will  be  subsequently 
referred  to  as  primitive  variables. 


The  second  set  is  (p  —  p7’u,/(7eA/g  Ag«un,  T^,  is  the  basic  state 

wall  temperature.  Assuming  that  the  equation  of  state  is  the  ideal  gets  law,  the  first 
disturbance  variable  is  proportional  to  the  disturbcince  temperature  and,  hence,  is 
zero  at  the  wall.  This  first  vauriable  is  used  instead  of  the  pressure  itself  because  the 
spectral  code  requires  that  the  boundary  conditions  for  each  vatriable  be  the  S3ime 
at  each  boundary.  With  this  restriction,  if  p  and  p  were  used  as  the  disturbance 
thermodynamic  variables  there  would  be  no  means  of  enforcing  the  wall  temperature 
boundary  condition. 

The  second  choice  of  dependent  variables  malces  the  inviscid  terms  in  the  origi¬ 
nal  governing  equations  homogeneous  functions  of  degree  one  and  the  viscous  terms 
homogeneous  functions  of  degree  zero.  (See,  for  instance,  Gelfand  and  Fomin  [1963] 
for  a  definition  of  homogeneous  functions  of  degree  k.  See  also  Schiff  and  Steger, 
[1979].)  These  variables  were  mainly  used  as  a  programming  stepping  stone  to  the 
third  choice.  They  will  be  subsequently  referred  to  as  conservative  variables. 

The  third  set  of  dependent  variables  is  y/gip  —  pT^I{'^eMl),p,  pv^,pv^,  pv^).  Here 
u^,  v^,  and  are  the  contravariant  velocity  components: 


dC 

=r  - U  + 

dx 


d(‘  dV  d(' 


(3.84) 


These  variables  will  be  subsequently  referred  to  as  scaled  conservative  variables. 

The  third  choice  is  available  made  because  the  continuity  equation  in  general 
curvilinear  coordinates  is  linear  in  these  variables.  In  fact,  all  of  the  coefficients  of  the 
dependent  variables  in  the  continuity  equation  are  constant.  Hence,  the  disturbance 
satisfies  conservation  of  mass  exactly  (and  uniformly)  even  when  the  parallel  flow 
approximation  is  made.  (The  parallel  flow  approximation  still  is  required  because  of 
nonparallel  terms  in  the  momentum  and  energy  equations.) 

The  main  input  required  for  each  of  these  sets  of  disturbance  variables  is  the  matrix 
Mj.  This  matrix  transforms  the  original  (default)  disturbance  dependent  variables  Q 
=  {Q^j  Q*,  Q^)  =  {p,  T,  u,  u,  iZ;)  into  the  new  dependent  veiriables: 


Q'  =  Afj(®,<)/^(^,  r)exp(e^(^^,^^,r))  complex  conjugate  (3.85) 


The  matrix  Mj  for  each  set  of  disturbance  variables  is  given  in  the  following 
subsections. 


Primitive  Variables 

In  this  case,  the  matrix  Mj  is  particularly  simple:  it  is  just  the  identity  matrix 


m; 


f  1  if  i  =  ;■ 
(  0  iii  ^  j 


All  of  its  spatial  and  temporal  derivatives  are  zero. 


(3.86) 
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Conservative  Variables 

The  matrix  Afj  for  the  conservative  variables  is 


/  1 


dp/ar 

0 

0 

0 


ikfkr) 

^  t£ 

P 

P 

_ ti; 

P 


0  0  0  \ 
0  0  0 
i  0  0 

p 

0  7  0 

oS  J 


(3.87) 


The  spatial  and  temporal  first  and  second  derivatives  of  this  matrix  can  be  computed 
using  the  chain  rule  and  the  derivatives  of  the  basic  state.  The  derivatives  of  the 
basic  state  aire  given  in  Section  3.2.3. 


Scaled  Conservative  Variables 
The  matrix  Afj  for  the  scaled  conservative  variables  is 

0  0  0  \ 

0  0  0 

\  dx  1  ^  I  dx 
pdi^  pdi^ 

1  ay  1  ay  1  dy 

p  p  p  04* 

1  az  I  at  10* 

^041  p04^  j 

(Note  that  even  if  the  grid  is  time  dependent  (i.e.,  dx/dr  ^  0,  etc.),  as  long  as 
r  ■=  t  and  there  is  no  disturbance  in  the  grid  speed,  the  disturbance  in  the  fourth 
contravau-iamt  velocity  component  =  1  -V  udr  fdx  vBt  jdy  xudr  fdz  is  identicaJly 
zero.)  The  spatial  and  temporail  first  and  second  derivatives  of  this  matrix  can  be 
computed  using  the  chadn  rule  and  the  derivatives  of  the  grid  and  baisic  state.  The 
derivatives  of  the  grid  were  given  eau’lier  in  Section  3.2.1.  The  derivatives  of  the  basic 
state  aire  given  in  Section  3.2.3. 


/  1 

apidp  1 

0/»/sr  0/»/0T 
0 

0 

0 


-( 


lezia)  ( 

9piaT)  V 


T... 

■ytM}. 


P 

~  P 
w 
P 


3.2.3  Basic  State  for  Parallel  and  Quasi- Parallel  Flows 

The  basic  state  and  its  spatial  and  temporal  derivatives  must  also  be  specified  for 
the  stability  analysis.  The  particular  vadues  assigned  depend  upon  whether  the  pair- 
allel  or  quasi-pairadlel  flow  approximation  is  made.  Also,  the  directions  in  which  the 
derivatives  are  talcen  depend  upon  the  coordinate  system  used.  Since  the  basic  state 
is  computed  using  the  similarity  coordinates,  its  derivatives  in  this  coordinate  system 
are  presented  first. 


24 


Basic  State  Derivatives  in  Similarity  Coordinates 

The  similarity  coordinate  system  has  already  been  described  in  Section  3.2.1.  Also, 
the  dependence  of  the  basic  state  on  these  coordinates  has  been  described  in  Sec¬ 
tion  3.1.1.  In  particular, 


p  = 

piv)Tiv)l{lM  =  l/(7eMn 

(3.89) 

T  = 

T{v) 

(3.90) 

P  = 

VTiv) 

(3.91) 

u  = 

m 

(3.92) 

V  = 

(3.93) 

w  = 

0 

(3.94) 

(See  equations  (3.13),  (3.25),  and  (3.26).)  Hence, 

(3.95) 

o 

11 

(3.96) 

-  =  0 
ae 

(3.97) 

dv  V 

^  “  ”2l 

(3.98) 

o 

11 

(3.99) 

and 

dp 

dr] 

0 

(3.100) 

11 

T'iv) 

(3.101) 

11 

fiv) 

(3.102) 

dv 

drj 

(3.103) 

dw 

dr) 

0 

(3.104) 

The  second  derivatives  are 

d'^p 

=  0 

(3.105) 
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d^T 

d^u 

W 

d^v 

d^w 

le 

d^P 

d^dp 

d^T 

d^dp 

d^u 

didp 

d^v 

d^ 

d^w 

didrf 

dp^ 

d^T 

drj^ 

d^u 

d^ 

572 


d^w 

dp^ 


0 

(3.106) 

0 

(3.107) 

Zv 

4^2 

(3.108) 

0 

(3.109) 

0 

(3.110) 

0 

(3.111) 

0 

(3.112) 

1  dv 

(3.113) 

2x  drj 

0 

(3.114) 

0 

(3.115) 

r(»/) 

(3.116) 

nv) 

(3.117) 

-  f'MT'M  -  fMr'M  ) 

(3.118) 

0 

(3.119) 

Basic  State  Derivatives  in  Rectilinear  Coordinates 

The  values  assigned  to  the  basic  state  and  its  derivatives  in  the  rectiline2ur  coordi¬ 
nate  system  depend  upon  whether  the  parallel  or  quasi-parallel  flow  approximation 
is  made. 

Parallel  Flow 
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When  the  basic  state  flow  is  approximated  as  being  locally  par2dlel,  it  depends  only 
upon  the  coordinate  j/;  it  is  independent  of  x.  Equivalently,  it  depends  only  upon  rj; 
it  is  independent  of  (See  equations  (3.63)  and  (3.64).)  Hence, 


Hence, 


P 

T 

u 

V 

w 


T{V) 

u{ti) 

0 


=  0 


dfj 

dfj 

du 

dfj 

dv 

dfj 

dw 

dfj 


[vVileMD  =  iK'leMl) 

(3.120) 

(3.121) 

(3.122) 

(3.123) 

(3.124) 

II 

o 

(3.125) 

o 

II 

1 

(3.126) 

o 

II 

3  I'Jt/- 

(3.127) 

^  =  0 
di 

(3.128) 

^  =  0 
d( 

(3.129) 

0 

(3.130) 

Ife).  r{rj) 

y/x/xoT{Tj) 

(3.131) 

Tiv)  . 

y/xlxoT{r]) 

(3.132) 

0 

(3.133) 

0 

(3.134) 

S  =  0 

de 

(3.135) 

d^T 

=  0 
de 

(3.136) 
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0 


(3.137) 


d^u 

de 


(3.138) 

(3.139) 


d^w 

dldfj 


(3.140) 

(3.141) 

(3.142) 

(3.143) 

(3.144) 


(3.145) 


(3.146) 


(3.147) 


(3.148) 
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0 


(3.149) 


d^w 

drj'^ 


Quasi-Parallel  Flow 

When  the  basic  state  flow  is  approximated  as  being  quasi-parallel,  all  of  the  coefl&- 
cients  in  the  disturbance  equations  which  are  functions  of  x  are  assumed  to  be  locally 
constant.  In  other  words,  they  are  retained  in  the  equations,  but  are  assumed  constant 
and  evaluated  at  the  particular  station  where  the  local  stability  analysis  is  performed. 


In  this  case, 

P  =  p(v)Tiv)/(l.Ml)  =  l/(leM^)  (3.150) 

T  =  T{ti)  (3.151) 

P  ==  l/Tiv)  (3.152) 

u  =  f{ri)  (3.153) 

u;  =  0  (3.155) 


(See  equations  (3.13),  (3.25),  and  (3.26).)  Applying  the  transformation  laws  (3.77), 
(3.78),  (3.79),  (3.81),  and  (3.83)  to  this  approximation  of  the  basic  state  yields: 


di 

du 

H 

dv 

dw 

dp 

dfj 

dfj 

du 

dfj 

dv 

dfj 


0 

(3.156) 

T\v) 

2xTiv)Jo 

(3.157) 

riv) 

2xT{ri)L 

(3.158) 

dv  1  r  mi  \J 

di  2xT{t})  Jo  ^  *  dp 

(3.159) 

0 

(3.160) 

0 

(3.161) 

(3.162) 

jy 

y/xlxoT{p) 

(3.163) 

W)  9v,^ 

^x/xoT{p)^n 

(3.164) 
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=  0 


(3.165) 


The  second  derivatives  are: 


(— r 

\2i7’(t;)  Jo 


de  ^  {2xT{Ti)f 


l-L-r 

\2xT{7])  Jo 


^  ^  ^  SS  Tjs^ds 
d'e  ~  de""  {2xT{Ti)f 


mri)  - 


(3.166) 


(3.167) 


2xT{ri) 


T(s)ds|  r\rf)  (3.167) 

/"'(,)  (3.168) 


(3.168) 


2  /•”  d^v 

Tiv)  Jo 


(3.170) 


d^p 

dldfj 

=  0 

(3.171) 

d^T 

dldfi 

m)  rn-?)  fn,)  p  i 

1  (3.172) 

T"(fi)  /-Vr.uJ 

2xT{7j)  Jo 

(3.173) 

d^u 

dldfj 

Tiv)  \nv)ir{v)  : 

[  (3.174) 

r  T(s)d3 

2xT{rj)  Jo 

(3.175) 
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dldff 


d^w 

didfj 


Tjrj)  \d^v  ^  I  i  T'ir,)  p, 


1  nrrf  s,  d'^v 

2xT{t))Jo 


d^p 

dfj^ 

d^T 

drj^ 


drj'^ 


=  0 


d^v 

dfi"^ 


d^w 

dfj'^ 
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\nfi) 

1  fmVriv) 

yJxlxQ 

^,\m)  r(,)J 

,  \  m  ] 

2 

^  rxv) 

\y/x^/XoT{Tj)} 

1 

1 

'T'ifj) 

1  (mynri) 

V^/xo 

T{v) 

rwj 

T'iv) 


f'M 


Tiv)  V 


y/xlxQT{rj)j 


nv) 


Inv) 

1 

frW)' 

l'  nr,)' 

dv 

T{v) 

yfxjxo 

\T{v)j 

T{ri) 

dp 

+ 


=  0 


[  T{rj)  Vd^v 
\y/x/xoT{Tj)f 


(3.176) 

(3.177) 

(3.178) 

(3.179) 

(3.180) 


(3.181) 

(3.182) 

(3.183) 


3.2.4  Computed  Results 

Some  results  are  presented  in  this  subsection.  They  were  obtained  using  the  Suther¬ 
land  viscosity  law  for  air,  constant  specific  heat  Cp  =  7e/(7e  —  1)  (with  7e  =  1.4), 
2uid  constant  Prandtl  number  =  0.70.  The  transformation  parameter  rjo  was  set  to 
4.5/^/2  3.181980516  in  order  to  make  comparisons  with  results  for  the  incompress¬ 

ible  boundary  layer.  ^ 


^  where  the  dehnition  of  t)  does  not  include  a  factor  of  y/^. 


Parallel  Flow,  Rectilinear  Coordinates 
Primitive  Variables 

A  special  insert  file  for  LSH  was  written  to  implement  the  equations  for  primitive 
variables  in  rectilinear  coordinates  in  order  to  make  direct  comparisons  between  its 
results  and  those  of  the  more  general  insert  file.  After  debugging,  both  of  them 
produced  identical  results  (since  they  use  the  same  formulation  of  the  equations). 

Particular  flows  investigated  are  some  of  those  published  in  the  literature  by  Malik 
(1990).  The  first  is  for  a  two-dimensional  disturbeince  (Malik’s  test  case  ifb):  M*  = 
10,  Te  =  lll.lllii^,  R  =  1000,  adiabatic  wall.  A  temporal  analysis  of  all  of  the 
eigenvalues  of  the  discretized  system  (using  31  collocation  points)  for  a  =  (0.12,0), 
=  (0, 0)  yields: 

uj  =  (0.1158601206,  1.367658245E-  04)  (3.184) 

Malik  reports  w  =  (0 . 1158647, 1.529E-04)  using  his  4CD  scheme  and  uj  = 
(0 . 1158519,1 .357E-04)  using  his  MDSP  method.  However,  Esfahtinieui  (1991),  who 
also  used  the  same  rectilinear  coordinates  for  his  stability  analysis  with  a  spectral  col¬ 
location  method  and  101  points,  reports  a  value; 

u  =  (0,115861436,  1.37657E-  04)  (3.185) 

which  agrees  much  better  with  our  results  as  well  as  with  Malik’s  MDSP  calculation. 
Having  obtained  w,  a  spatial  eigenvalue  calculation  yields: 

w  =  (0.1158,0) ,  Q  =  (0.11994036,  -1.3877400E-  04)  (3.186) 

along  with  the  corresponding  eigenvector.  The  PSE  code  PSX  reproduces  this  result 
when  mairched  in  a  parallel  mode  with  this  as  the  initial  condition.  A  comparison^ 
of  the  derivative  of  the  linear  operator  used  in  the  locad  anedysis  with  respect  to  or 
and  the  coefficient  of  the  derivative  of  the  shape  function  with  respect  to  x  in  the 
PSX  code  also  yields  complete  agreement.  (The  analytical  derivative  used  in  the  local 
analysis  also  agrees  very  well  with  second  order  accurate  finite  differences.) 

Results  for  Malik’s  test  case  #3  for  a  three-dimensional  disturbance  have  also 
been  obtained.  This  test  case  is  for  an  adiabatic  wall  with  a  boundary  layer  edge 
Mach  number  Me  =  2.5.  The  edge  temperature  is  Tg  =  148.148148iif.  The  Reynolds 
number  is  R  =  3000.  Again,  31  collocation  points  were  used.  A  temporal  analysis  for 
the  entire  spectrum  of  the  discretized  system  for  a  =  (0.06,0),  ^  =  (0, 10,0)  yields 
as  one  of  the  eigenvalues: 

u  =  (3.673723764E  -  02,  5.895818712E  -  04)  (3.187) 

Malik  reports  uj  =  (3.67340E  —  02, 5.840E  —  04).  LSH  gives  as  the  corresponding 
spatial  mode  for  u)  =  (3.67E  —  02, 0),  0  =  (0.10,0): 

a  =  (5.9977439E  -  02,  -8.0711300E  -  04)  (3.188) 

^  using  the  UNIX  command 
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Mau^ching  with  the  PSE  code  in  a  parallel  flow  mode  using  the  spatial  eigenvalue  and 
eigenvector  ais  initial  conditions  reproduces  the  same  result.  Also,  the  coefficients 
of  the  derivatives  of  the  local  lineau:  operator  with  respect  to  q  and  d  in  LSH  again 
agree  very  well  with  the  derivatives  of  the  shape  function  with  respect  to  i  and  2, 
respectively,  in  the  PSE  code  (as  well  as  finite  difference  calculations). 

Our  final  test  case  is  a  comparison  with  the  local  stability  results  produced  by  an 
independently  developed  fourth  order  accurate  compact  finite  difference  code  written 
at  DynaFlow.  This  code  uses  the  Euler-MacLaurin  finite  difference  formula  (e.g., 
Malik  1990). 

The  conditions  are  the  same  as  those  studied  by  El-Hady  (1980)  a’-d  Bertolotti 
(1991).  The  boundary  layer  edge  Mach  number  is  Me  =  1-6.  The  edge  emperature 
is  Te  =  206 A",  eind  the  Reynolds  number  is  R  =  400. 

When  uj  =  (0.016,0),  3  =  (0,0),  and  the  outer  computational  domain  is  tidcen  to 
be  T}rnar  =  100,  the  compact  finite  difference  code  yields: 

a  =  (3.509E  -  02,  8.63E  -  04)  (3.189) 

When  Tjmax  =  200,  it  yields: 

a  =  (3.5E-02,  8.76E-04)  (3.190) 

Using  51  collocation  points  to  compute  the  same  spatial  eigenvalue,  LSH  gives 

a  =  (3.5187357E  -  02,  8.7718022E  -  04)  (3.191) 

The  accuracy  of  the  spatial  growth  rate  substantially  deteriorates  if  only  31  collocation 
points  are  used.  This  Ccise  is  considered  further  in  the  following  section. 

Conservative  Variables* 

The  last  test  case  described  above  has  also  been  used  to  check  the  correctness  of 
the  insert  files  when  the  conservative  variables  are  used.  Again,  this  is  a  Mach  1.6 
boundary  layer  edge  flow  with  an  edge  temperature  of  206A’.  u;  =  (0.016, 0)  and 
3  =  (0,0).  Since  the  discretization  error  is  different  than  for  the  primitive  variables, 
the  grid  can  be  consecutively  refined  and  the  results  compared  to  see  whether  they 
approach  one  another: 


no.  pts. 

Q  =  (ar,Q,) 

primitive 

conservative 

(3.5187357E-02,8.7718022E-04) 

(3.5190580E-02,8.767l680E-04) 

(3.5190440E-02,8.7695450E-04) 

(3.5190019E-02,8.7687392E-04) 

(3.5187377E-02,8.7707594E-04) 

(3.5190577E-02,8.7673937E-04) 

(3.5190431E-02,8.7696982E-04) 

(3.5190015E-02,8.7687837E-04) 

The  results  seem  to  be  converging,  albeit  slower  them  one  might  expect.  We 
conclude  that  the  equations  insert  files  pass  this  test. 

*One  must  carefully  specify  the  boundeuy  conditions  in  the  Definitions  file  in  this  case.  Here, 
the  first  disturbance  variable  must  be  zero  at  the  wall.  Also,  there  is  no  boundary  condition  on  the 
second  disturbance  variable  (p). 
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Quasi-Parallel  Flow,  Rectilinear  Coordinates 
Primitive  Variables 

Continuing  with  the  same  test  case  described  above,  the  quasi-parallel  flow  approxi¬ 
mation  yields: 

a  =  (3.4436138E  -  02,  1.0105549E  -  03)  (3.192) 

using  121  collocation  points. 

Conservative  Variables 

Again,  for  the  same  boundary  layer  (adiabatic  wall,  Mg  =  1.6,  Te  =  206K),  the 
quasi-parallel  flow  approximation  yields: 

a  =  (3.4660S56E  -  02,  1.1219040E  -  03)  (3.193) 

using  the  conservative  variables  with  ui  =(0.016,0),  /3  =(0,0)  and  121  collocation 
points.  This  obviously  is  in  slight  disagreement  with  (3.192).  This  is  precisely  because 
the  basic  state  depends  on  x.  In  particular,  the  transformation  between  the  primitive 
and  conservative  variables  depends  upon  u  and  p,  which  in  turn  depend  upon  x  for 
any  fixed  value  of  y.  Hence,  it  is  not  possible  for  both  sets  of  variables  to  be  normad 
modes,  yet  this  is  what  is  assumed  by  the  local  analysis.  The  disturbance  equations 
are,  thus,  slightly  different,  leading  to  slightly  different  results.  Additional  support 
for  this  can  be  found  in  the  following  subsection. 

Quasi-Parallel  Flow,  Similarity  Coordinates 
Primitive  Variables 

Again,  continuing  with  the  adiabatic  wall  Mach  1.6  boundary  layer  (with  edge  tem¬ 
perature  Te  =  206X),  one  obtains 

a  =  (3.4967897E  -  02,  9.0065327E  -  04)  (3.194) 

for  the  two-dimensional  disturbance  w  =  (0.016,0),  =  (0,0)  using  51  collocation 

points,  and 

a  =  (3.4969567E  -  02,  8.9995719E  -  04)  (3.195) 

using  121  collocation  points.  This  last  result  should  be  compared  to  that  obtained 
using  the  conservative  variables. 

Conservative  Variables 

For  the  same  test  case,  using  the  conservative  variables  and  121  collocation  points 
gives: 

a  =  (3.4967547E  -  02,  8.9993425E  -  04)  (3.196) 

This  is  virtually  the  same  as  (3.195).  This  can  be  explained  as  follows.  In  the  similar¬ 
ity  coordinates,  u  and  p  depend  upon  rj  only.  The  normal  velocity  component  (which 
is  very  small)  does,  however,  depend  weakly  upon  x.  Hence,  the  transformation  be¬ 
tween  the  primitive  and  conservative  v^iables  is  (almost)  independent  of  x.  This,  in 
turn,  means  that  the  variables  should  have  (Adrtually)  the  same  spectrum. 
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Scaled  Conservative  Variables 

Continuing  with  the  same  test  case,  but  using  51  collocation  points  for  greater  effi¬ 
ciency,  one  obtains 

a  =  (3.4966572E  -  02,  -1.5638850E  -  04)  (3.197) 


using  the  scaled  conservative  variables.  This  is  approximately  what  one  should  expect 
because,  from  (3.47), 

Hence, 


di 


(3.4966572E-  02,  1.09361150E-  03) 


(3.199) 

(3.200) 

(3.201) 

(3.202) 

(3.203) 


This  compares  better  with  (3.194). 

Again,  the  spatial  eigenvalue  (3.197)  can  be  reproduced  by  the  PSE  code  by  march¬ 
ing  in  a  parallel  mode  starting  with  the  correct  eigenvalue  and  eigenvector.  Also,  the 
UNIX  debugger  has  been  used  to  check  that  the  coefficients  of  the  disturbance  quan¬ 
tities  in  the  conservation  of  mass  equation  are  correct  when  the  scaled  disturbance 
variables  are  used.  This  is  a  signific2uit  (although  not  comprehensive)  check  because 
of  the  many  intermediate  algebraic  manipulations  that  the  code  performs  to  achieve 
the  final  (in  this  case,  simple)  result. 


3.2.5  Linear  PSE 

The  same  flat  plate  boundary  layer  considered  in  the  previous  sections  was  also  an¬ 
alyzed  with  the  linear  PSE.  Again,  this  is  a  Mach  1.6  boundary  layer  edge  flow  with 
an  edge  temperature  of  206K,  u  =  (0.016,0)  and  0  =  (0,0).  The  results  obt2dned 
with  at  step  size  =  30  are  shown  in  Figure  3.1.  Similar  results  are  obtained  when 
a  step  size  of  20  is  used,  but  with  a  step  size  of  10  the  computed  growth  rate  begins 
to  oscillate  in  the  streamwise  direction  due  to  the  pressure  gradient  dp/d^. 

The  linear  PSE  integration  used  approximate  initial  conditions  obtained  from 
a  quasi-par2iliel  local  analysis  in  similarity  coordinates.  Figure  3.1  compares  these 
results  with  those  obtained  from  another  program  which  uses  the  compact  Euler- 
MacLaurin  finite  difference  formula  (i.e.,  Malik,  [1990])  to  solve  the  PSE  specifically 
for  the  flat  plate  compressible  boundary  layer.  The  transients  for  the  two  PSE  calcu¬ 
lations  are  different  because  the  initial  conditions  were  different:  the  compact  finite 
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difference  program  used  initial  conditions  generated  from  a  local  quasi-p^lraIlel  anal¬ 
ysis  in  Cartesian  coordinates.  The  amplitudes  shown  were,  thus,  scaled  so  that  the 
last  data  point  computed  by  the  compact  scheme  agreed  with  the  spectral  calculation. 
This  maJces  it  clearer  that  the  differences  ate  due  to  the  transients. 

3.2.6  Nonlinear  PSE 

Finally,  the  nonlinear  PSE  were  used  to  predict  the  spatied  evolution  of  a  subharmonic 
disturbance  in  a  laminar  flow  modulated  by  a  two-dimensionaJ  primary  instability. 
The  boundary  conditions  are  the  same  as  those  of  El-Hady  (1991),  with  whom  we 
compare  our  results.  The  edge  temperature  is  300/^,  the  Sutherlemd  viscosity  law  is 
used,  the  Prandtl  number  is  0.70,  and  the  edge  Mach  number  varies  from  0  to  1.2.^ 
In  each  case,  the  mean  flow  distortion  and  hiirmonics  created  by  the  primary  wave 
are  neglected.  The  primary  wave  itself  is  assumed  to  grow  lineaurly. 

The  results  for  an  edge  Mach  number  =  1.2  are  shown  in  Figures  3.2  zmd  3.3. 
The  frequency  of  the  primary  disturbance  is  F  =  IO^lj/R  is  60  emd  the  spanwise 
wavenumber  of  the  subharmonic  is  6  =  IO^0/R  is  0.15.  R  is  the  square  root  of  the  x 
Reynolds  number. 

Two  different  PSE  results  are  shown  in  Figures  3.2  amd  3.3.  The  solid  lines 
illustrate  the  results  obtained  when  Cartesian  coordinates  aure  used,  amd  the  daished 
lines  show  the  results  obtained  when  similarity  coordinates  are  used.  The  x  and  * 
symbols  represent  El-Hady’s  results. 

Figure  3.2  compares  the  spatiad  growth  rate  predictions  at  the  local  maiximum  in 
the  rms  profiles.  One  can  see  that,  apart  from  some  initial  tramsients,  the  results 
using  Cartesian  and  similarity  coordinates  agree  very  well.  This  is  an  importamt  test 
because  the  Cartesian  coordinates  are  orthogonal  whereas  the  similairity  coordinates 
are  not.  The  agreement  with  El-Hady’s  results  is  not  ais  good,  however.  This  is 
commented  on  in  further  detail  later. 

Figure  3.3  compares  the  values  of  ^{89 /d(^)  =  a,  predicted  by  the  PSE  in  Carte¬ 
sian  and  similarity  coordinates  to  El-Hady’s  results.  Since  — q,  is  readly  the  growth 
rate  of  the  norm  of  the  variable  which  is  used  to  partition  the  dependence,^  amd  all 
of  the  disturbances  varnish  in  the  far-field,  it  cam  be  shown  that  the  norm  should  be 
the  same  regardless  of  whether  Cartesiam  or  similau-ity  coordinates  aure  used.  While 
the  agreement  between  the  Cartesian  amd  similarity  results  in  Figure  3.3  is  not  as 
good  ais  the  local  results  shown  in  Figure  3.2,  it  still  appeal’s  reasonable  given  the  size 
of  the  transients. 

Some  of  the  discrepamcies  between  the  PSE  results  amd  El-Hady’s  results  may  be 
due  to  the  following.  While  the  PSE  calculations  neglect  the  derivatives  and  the 
normal  velocity  component  of  the  basic  state  to  better  simulate  El-Hawly’s  theoretical 
results,  they  retain  the  derivatives  of  the  mode  shape  functions.  This  introduces  a 
nonlocal  effect  which  is  not  present  in  El-Hady’s  results.  In  fact,  El-Hady’s  analysis 
assumes  that  the  basic  flow  modulated  by  the  primajy  disturbamce  is  locally  a  periodic 
function  of  time  r  and  distance  Hence,  it  neglects  the  spatial  growth  rate  of  the 

^  El-Hady  also  investigated  some  other  edge  Mach  numbers. 

^here  taken  to  be 
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Figure  3.2:  Spatial  Growth  Rate  at  Maximum  in  rms  Profiles.  Amplitude  of  primary 
disturbance  at  neutral  point  =  0.1%.  Frequency  F  =  60,  spanwise  wavenumber 
h  =  0.15.  □:  PSE  primary  mode  (2,0),  O:  PSE  subharmonic  mode  (1,1).  Solid 
line:  Cartesian  coordinates.  Dashed  line:  similarity  coordinates,  x:  El-Hady  (1991) 
primary  mode  (2,0),  *:  El-Hady  (1991)  subharmonic  mode  (1,1). 

primary  disturbance.  The  PSE  calculations,  however,  include  the  spatial  growth  rate 
of  the  primeu'y  disturbance.  The  transport  and  thermodynamic  properties  are  also  dif¬ 
ferent  than  what  El-Hady  used.  El-Hady  computed  the  viscosity  from  the  Sutherland 
law.  However,  he  determined  the  thermal  conductivity  not  from  the  constant  Prandtl 
number  and  specific  heat  but  from  k*  =  0.6325x10“®  V^/(l-|-(245.4/r)xl0“^*),  where 
k  is  meaisured  in  calories  per  centimeter  per  second  per  degrees  Celsius  2«id  T  is  in 
Kelvin.  He  also  used  the  National  Bureau  of  Standards  perfect-g2is  tables  to  spec¬ 
ify  the  Prandtl  number.  He  then  computed  the  specific  heat  and  enthailpy  from  the 
Prandtl  number  knowing  the  viscosity  and  thermal  conductivity. 


3.3  Summary 

After  having  performed  the  following  tests: 

1.  Manually  cross-checked  the  calculation  of  the  metric  terms. 

2.  Computed  some  of  the  differential  operators  two  different  ways. 

3.  Compared  with  results  available  in  the  literature. 

4.  Compared  with  calculations  of  different  software  written  at  DynaFlow  by  dif¬ 
ferent  personnel. 
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Figure  3.3:  Spatial  Growth  Rate  at  Maximum  in  rms  Profiles.  Amplitude  of  primary 
disturbance  at  neutral  point  =  0.1%.  Frequency  F  =  60,  spanwise  wavenumber 
h  =  0.15.  □:  PSE  Primary  mode  (2,0),  O:  PSE  Subharmonic  mode  (1,1).  Solid 
line:  Cartesian  coordinates.  Dashed  line:  similarity  coordinates,  x:  El-Hady  (1991) 
primary  mode  (2,0),  *:  El-Hady  (1991)  subharmonic  mode  (1,1). 

5.  Compared  PSE  results  obtained  in  Cartesian  (orthogon^ll)  and  similarity  (non- 
orthogoncil)  coordinates. 

the  suite  of  local  linear  stability  analysis  aind  PSE  software  appe^u^s  to  be  working  as 
designed.  This  should  help  give  more  confidence  in  the  following  applications  which 
have  more  relevance  to  hypersonic  flight. 
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Chapter  4 

Hypersonic  Flow  over  a  Sharp 
Cone  at  Angle  of  Attack 


The  stability  of  the  Mach  8  flow  past  a  sharp  7*  half-angle  cone  at  angle  of  attack  has 
been  studied  experimentally  by  Stetson  et  al.  (1985,1992).  They  measured  the  mass 
flux  and  total  temperature  fluctuations  along  the  windward  and  leeward  meridiams  at 
2°  angle  of  attack  but  only  along  the  windward  meridian  at  4°  angle  of  attack.  The 
free-stream  static  temperature  was  approximately  95®  i?,  and  the  surface  temperature 
was  allowed  to  reach  equilibrium  (adiabatic  wall).  Finally,  the  free-stream  Reynolds 
number  of  the  flow  for  the  windward  meridian  measurements  was  1  million  per  foot. 
However,  they  lowered  the  free-stream  Reynolds  number  to  0.5  million  per  foot  for 
their  measurements  along  the  leeward  meridian  to  obtain  a  greater  length  of  laminar 
flow. 

The  nximerical  analyses  presented  here  simulate  the  conditions  of  their  experiment 
at  2®  angle  of  attack.  However,  the  unit  Reynolds  number  is  0.5  million  per  foot  for 
both  the  windward  and  leeward  meridians  and,  hence,  different  than  the  value  for  the 
windward  meridian  experiment.  Since  the  unit  Reynolds  number  has  no  beauring  on 
the  present  theoretical  results,  this  will  not  affect  the  comparison  of  the  dimensionless 
results  with  the  experimental  data. 

Only  limited  linear  local  and  nonlocal  numerical  results  were  obtained  because  of 
the  unexpected  size  of  the  nonparallel  effects.  These  nonparallel  effects  are  so  large 
that  we  spent  a  good  deal  of  time  confirming  them  by  performing  a  number  of  unique 
tests  which  are  documented  below.  The  tests  reveal  that  the  stability  codes  produce 
self-consistent  results  and  that  the  nonp2krallel  effects  are  indeed  large. 


4.1  Mean  Flow  Calculation 

The  mean  flow  was  calculated  using  the  AFWAL  Parabolized  Navier-Stokes  (PNS) 
code.  This  code  h£is  been  developed  by  a  number  of  different  workers  since  1979  and 
has  been  documented  by  Kaul  and  Chaussee  (1983);  Stalnaker,  Nicholson,  Hanline, 
and  McGraw  (1986);  and  Rajendran  (1989).  It  uses  a  space-marching  technique  and 
body-fitted  coordinates  to  integrate  the  PNS  equations  downstream  of  data  specified 
in  an  initial  axial  plane  =  Xq. 
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4.1.1  Coordinate  System 

The  permissible  body-fitted  coordinates  are  all  of  the  form: 


=  x\l)  (4.1) 

=  x\j,K)  (4.2) 

=  x%J,K)  (4.3) 

where  is  the  body  axis.  (See  Figure  4.1.)  This  means  that  the  solution  is  marched 


Figure  4.1:  Coordinate  System  Used  in  the  Stability  Analyses. 

in  planes  perpendicular  to  the  body  axis.  Hence,  the  body-fitted  coordinates  are 
nonorthogonal  whenever  the  slope  of  the  body  in  the  axial  direction  is  nonzero.  This 
is  certaunly  the  case  for  any  cone  of  nonzero  half-angle.  We  note  this  here  because 
virtually  the  same  coordinate  system  (Section  4.2)  was  used  for  the  stability  analyses 
in  order  to  simplify  the  problem  of  interpolating  the  basic  state  solution  onto  the 
disturbance  grid. 

4.1.2  Initial  Conditions 

Two  methods,  each  applicable  to  a  different  type  of  geometry,  can  be  used  to  generate 
initial  data  for  the  PNS  equations.  One  of  these  is  the  “step-back”  procedure  of  Schiff 
and  Steger  (1979).  Another  is  the  solution  of  the  governing  equations  in  the  forebody 
region. 

The  “step- back”  procedure  of  Schiff  and  Steger  (1979)  has  been  implemented  in 
the  AFWAL  PNS  code.  It  can  be  used  for  pointed  forebodies  which  have  an  almost 
conical  flowfield  about  them.  In  this  case,  a  conical  coordinate  system  is  chosen 
which  includes  the  initial  plane  and  a  guess  is  made  for  the  initial  data.  A  step  is 
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then  taken  in  the  marching  direction  to  obtain  the  solution  immediately  downstream. 
Approximating  the  solution  ais  locally  conical,  the  solution  in  the  downstream  plane 
should  be  the  saune  as  that  in  the  initial  plane  along  the  same  conicai  rays.  Hence, 
the  initial  plane  of  data  is  overwritten  with  the  solution  downstream  adong  the  same 
conical  rays.  This  yields  a  new  guess  for  the  initial  data  which,  upon  stepping  back  to 
the  initial  station,  can  be  used  to  begin  the  integration  anew.  Repeating  this  process 
until  the  gradients  of  the  solution  vanish  along  the  conical  rays  to  some  prescribed 
level  gives  the  locally  conical  flow.  This  is  how  the  initial  data  were  generated  for  the 
sharp  cone  studies  presented  here. 

For  blunt  bodies,  a  different  solution  methodology  must  be  used  because  the 
governing  equations  in  the  nose-tip  region  are  elliptic.  A  thin-layer  Navier-Stokes 
(TLNS)  code  developed  by  Rizk,  Scott,  and  Newman  (1983)  can  and  has  been  used 
for  this  purpose.  In  particular,  their  code  was  used  to  generate  the  starting  data 
for  the  blunt  cone  PNS  solution  used  in  Section  5.  However,  this  TLNS  code  uses 
a  variation  of  spherical  coordinates  while  the  AFWAL  PNS  code  uses  a  variation 
of  cylindrical  coordinates.  Multidimensional  interpolation  must,  thus,  be  done  to 
transfer  the  data  from  the  TLNS  grid  to  the  AFWAL  PNS  initial  plane.  Routines 
supplied  with  the  AFWAL  PNS  code  were  used  for  this  purpose. 

4.1.3  Basic  State  Results 

In  all  the  PNS  calculations,  the  equations  were  first  transformed  to  cylindrical  coordi¬ 
nates  and  then  to  computational  coordinates  I,J,K.  In  the  computational  coordinates, 
the  cone  surface  and  bow  shock  appear  as  the  surfaces  K  =  I  and  K  =  KMAX^ 
respectively.  The  251  points  were  used  in  the  radial  direction,  K,  and  19  points  in  the 
circumferential  direction,  J.  The  equations  were  then  discretized  amd  solved  with  the 
Beam- Warming  algorithm  in  conjunction  with  the  Vigneron  technique  to  suppress  the 
“departure  solutions.”  The  gas  was  assumed  perfect  and  the  cone  adiabatic.  Finally, 
PL0T3D  FORTRAN  unformatted  planes  files  (Walatka,  et  al.,  1990)  were  created  and 
used  to  inspect  the  computed  solution.  These  PL0T3D  files  were  then  input  to  the 
stability  codes  for  the  stability  analyses. 

Some  of  the  basic  state  windward  profiles  and  their  first  and  second  derivatives 
with  respect  to  eire  shown  in  Figures  4.2  through  4.4. 

The  length  scale  for  is  y/Rei^  where  Rcl  is  the  Reynolds  number  based  on 
the  free-stream  shock  conditions  and  length  L  =  40”  of  the  cone.  The  scales  for  the 
dependent  variables  are  all  taken  to  be  the  corresponding  free-stream  values.  The 
bow  shock  is  at  =  26.7  in  each  figure,  very  near  the  last  plotted  point. 

The  profiles  have  been  obtained  from  the  PNS  data  by  interpolating  them  with 
fifth  order  polynomials  in  and  Most  of  the  data  are  reasonably  smooth, 

but  the  same  anomalies  seen  in  the  blunt  cone  contravariant  velocity  profiles  near 
the  wall  (Figure  4.3  b)  are  also  seen  here. 

The  leeward  interpolated  profiles  at  R  =  712  are  shown  in  Figures  4.5  through 
4.7. 

The  bow  shock  is  at  —  37.44  at  this  station,  much  further  from  the  body  than 
it  is  on  the  windward  meridian. 
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Note  that  the  boundaxy  layer  is  approximately  twice  <is  thick  on  the  leeward 
meridiem  as  it  is  at  the  seime  station  on  the  windward  meridian.  This  is  because  the 
boundaxy  layer  is  compressed  much  more  on  the  windward  side  by  the  higher  pressure 
behind  the  stronger  bow  shock  created  by  the  greater  local  deflection  of  the  flow.  In 
addition  to  the  boundaxy  layer  being  thicker  on  the  leewaxd  side,  the  corresponding 
basic  state  derivatives  axe  also,  in  general,  smaller.  Finally,  tne  computed  derivatives 
appear  to  have  some  small  numerical  oscillations  in  them  in  the  lower  hadf  of  the 
boundary  layer.  However,  no  effort  was  made  to  improve  upon  these  results  by 
rerunning  the  AFWAL  PNS  code  as  our  experience  with  the  blunt  body  calculation 
indicated  that  this  was  very  difficult  to  do. 


4.2  Disturbance  Coordinate  System 

The  disturbance  coordinate  system  is  nonorthogonal  and  almost  identical  to  the  al 
gebrciic  one  that  the  AFWAL  PNS  code  uses.  It  is  given  by; 

x'^  =  i2sin(^^  — tt) 

=  Rcos{^^  —  If) 

R  =  n{i\e)+e 

where  is  the  distance  measured  along  the  body  axis,  is  the  circumferential  angle 
measured  in  radians  from  the  windward  meridian  in  the  — direction,  rj,  is  the  local 
body  radius,  and  is  the  radial  distance  from  the  body  surface.  See  Figure  4.1.  The 
similarity  between  the  AFWAL  PNS  algebraic  grid  and  the  coordinate  system  used 
for  the  stability  analysis  is  intentional  and  enables  the  basic  state  to  be  interpolated 
with  greater  ease  and  accuracy. 


-J?sin(^^) 

-i?cos(^^) 


(4.4) 


4.3  Local  Linear  Stability  Analyses 

As  mentioned  earlier,  the  nonparallel  effects  for  the  sharp  cone  at  angle  of  attack 
axe  so  large  (see  Section  4.4)  that  we  performed  a  number  of  tests  to  check  the 
consistency  of  the  results.  Virtually  all  of  these  tests  were  conducted  for  the  flow 
along  the  windward  meridian. 

4.3.1  Windward  Meridian 

One  consistency  check  we  performed  was  that  the  computed  eigenvalue  should  be  the 
same  regardless  of  whether  the  velocity  components  or  the  mass  fluxes  axe  used  as 
dependent  variables  in  the  disturbance  equations  as  long  as  all  of  the  derivatives  of 
the  basic  state  flow  with  respect  to  and  are  neglected}  When  these  derivatives 
are  neglected,  the  density  is  approximated  as  being  independent  of  and  and  so 
the  transformation  matrix  between  the  different  dependent  vciriables  is  independent  of 

^and  as  long  as  the  truncation  errors  in  the  solution  of  the  disturbance  equations  are  negligible. 
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and  In  this  czise,  the  assumption  of  norm^d  modes  for  the  velocity  components 
is  equivalent  to  the  assumption  oi  normal  modes  for  the  mass  fluxes.  However,  if  the 
transformation  matrix  is  a  function  of  any  independent  variable  in  which  the  stability- 
equations  axe  assumed  to  be  separable,  the  normal  mode  approximation  for  the  pre- 
and  post-transformed  dependent  variables  is  different}  In  this  case  an  inconsistency 
appears  in  the  equations  which  cannot  be  avoided. 

That  the  eigenvalues  axe  the  same  when  the  and  derivatives  of  the  b«isic 
state  flow  are  neglected  has  been  confirmed.  Some  supporting  results  are  shown  in 
Figure  4.8.  Figure  4.8  shows  the  local  spatial  growth  rate  variation  of  a  70kHz  second 
mode  \  ‘  isus  the  local  Reynolds  number  b«ised  on  x^Jy/R^.  Although  there  appear 
to  be  only  three  curves  in  Figure  4.8,  there  axe  in  fact  four.  Two  of  the  curves  axe 
•irtuaJly  coincident.  The  coincident  curves  are  the  results  obtained  from  the  described 
test. 

While  the  results  shown  agree,  this  test  revealed  that  the  interpolation  must  be 
done  very  carefully.  When  the  basic  state  derivatives  were  determined  by  differ¬ 
entiating  the  interpolating  polynomials,  the  eigenvadues  agreed.  Other  interpolation 
schemes  that  used  central-differenced  derivatives  at  the  basic  state  grid  points  and 
more  compact  higher  order  interpolants  gave  eigenvalues  differing  by  about  10%. 

When  the  and  derivatives  of  the  basic  state  are  included,  they  lead  to  quasi¬ 
parallel  tern,s  in  the  equations  which  are  different  in  each  dependent  variable  formu¬ 
lation.  The  computed  solutions  are  then  different.  The  two  quasi- paurallel  solutions 
obtained  from  these  different  formulations  axe  also  shown  in  Figure  4.8.  The  discrep¬ 
ancy  in  the  quasi-parallel  results  is,  thus,  some  indication  of  the  inaccuracy  of  the 
parallel  flow  approximation. 

Finally,  although  all  of  the  results  shown  in  Figure  4.8  vary  smoothly  with  R, 
they  may  not  if  the  interpolation  is  not  of  high  enough  order.  The  results  shown 
in  Figure  4.8  were  all  obtained  by  interpolating  the  basic  state  solution  with  fifth 
order  Lagrange  polynomials  in  and  However,  if  the  order  of  the  Lagrange 

interpolants  in  and  is  lowered  to  1  (linear),  some  of  the  results  are  not  as 
accurate.  Refer  to  Figure  4.9. 

Figure  4.9  shows  the  results  obtained  when  the  derivatives  of  the  baisic  state  with 
respect  to  and  are  neglected.  Here,  the  order  of  the  £md  interpolating 
polynomials  does  not  seem  to  have  much  of  an  effect  on  the  accuracy.  However,  when 
the  derivatives  of  the  basic  state  with  respect  to  and  are  included  (Figure  4.10), 
the  higher  order  interpolants  are  needed  to  obtain  accurate  results.  This  is  especially 
true  if  the  order  of  the  Lagrange  interpolant  in  is  reduced  to  3  from  5.  See 
Figure  4.11.  The  higher  order  accuracy  is  necessitated  by  the  large  effect  that  the 
derivatives  of  the  basic  state  have  on  the  stability  analysis  and  which,  thus,  must  be 
computed  accurately.  This  requirement  also  took  some  time  to  identify. 

The  tests  described  above  indicate  that  self-consistent  results  may  be  obtained 
with  the  local  stability  code  ais  long  as  consistent  approximations  are  made.  This 

^Another  example  is:  the  normal  mode  approximation  for  the  Cartesian  velocity  components  is 
not  the  same  as  the  normal  mode  approximation  for  the  curvilinear  velocity  components  because 
the  direction  and  magnitude  of  the  basis  vectors  of  the  curvilinear  coordinate  system  are  in  genera, 
functions  of  position. 
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holds  true  not  only  for  the  formulation  of  the  problem  via  the  selection  of  a  set 
of  dependent  variables,  but  also  for  the  basic  state  interpolation  scheme.  That  the 
results  agree  in  the  final  analysis  helps  validate  our  work. 

Similar  validation  tests  of  the  PSE  software  were  also  performed  and  axe  presented 
in  Section  4.4.  Before  proceeding  to  them,  however,  we  present  the  locad  linear  results 
we  obtciined  on  the  leeward  side  of  the  cone. 

4.3.2  Leeward  Meridian 

More  limited  results  were  obtained  on  the  leeward  side  as  the  contract  period  was 
ending.  Figure  4.12  shows  the  spatial  growth  rates  obtained  from  different  local 
amalyses.  Two  of  the  local  analyses  include  curvature  but  exclude  the  derivatives  of 
the  basic  state  with  respect  to  and  They  use  different  dependent  variables, 
though.  The  remaining  two  analyses  include  the  basic  state  derivatives  in  a  quasi¬ 
parallel  manner. 

As  expected,  the  results  obtained  are  the  same  regardless  of  which  dependent 
variables  are  used  as  long  as  the  basic  state  derivatives  with  respect  to  and 
are  neglected.  However,  the  results  are  very  different  when  these  derivatives  are 
included  in  a  quasi-parallel  manner.  Moreover,  although  not  shown,  the  predicted 
values  of  the  real  part  of  Or  also  differ  much  more  than  usual  (about  5%).  These 
large  differences  make  it  difficult  to  determine  whether  the  results  all  correspond  to 
the  same  family  of  disturbances.  Preliminary  results  seem  to  indicate  that  there  may 
be  other  neighboring  families  of  disturbances. 


4.4  Linecir  PSE  Analyses 

As  mentioned  earlier,  the  nonparallel  effects  that  we  have  found  are  unusually  large  in 
the  flows  along  the  windward  <ind  leeward  meridians.  Similar  testing  of  the  windward 
meridian  results  with  different  sets  of  dependent  variables  was  done  to  confirm  that 
these  nonparallel  effects  are  real.  In  each  test,  the  partition  (2.47)  between  the  shape 
function  and  exponential  variation  was  applied  to  the  same  variable:  the  disturbance 
density. 

4.4.1  Windward  Meridian 

As  with  the  local  stability  analyses,  we  checked  our  results  using  two  different  sets  of 
variables:  f,  T,  and  the  physical  contravariant  velocity  components]  and  p,  p,  and  the 
physical  contravariant  mass  fluxes.  However,  no  longer  is  the  disturbance  assumed 
to  depend  solely  on  the  local  properties  of  the  flow  profiles.  The  modulations  of 
the  mode  amplitude  and  phase  due  to  the  nonparaillel  basic  flow  are  now  talcen  into 
account. 

If  the  basic  state  derivatives  with  respect  to  and  are  neglected  while  the  true 
(nonparallel)  profiles  are  used  in  integrating  the  disturbance  downstream,  the  results 
shown  as  dotted  lines  in  Figure  4.13  are  obtained.  Although  the  results  obtained 
using  different  dependent  variables  are  initially  the  same  as  discussed  in  Section  4.3, 
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they  diverge  downstream.  Note  that  this  would  not  be  the  case  if  the  basic  state  flow 
were  actually  parallel  for  then  both  results  would  be  equivalent  to  the  local  stability 
results  which  are  identical.  Hence,  it  appears  that  the  derivatives  of  the  basic  state 
with  respect  to  and  are,  in  fact,  not  necesstirily  negligible. 

If  the  profiles  axe  allowed  to  vary  and  the  first  derivatives  of  the  basic  state  with 
respect  to  and  are  retained  in  the  equations,  then  the  differences  vanish.  This 
is  shown  by  the  agreement  of  the  two  solid  lines  in  Figure  4.13  further  downstream. 
These  lines  initially  disagree  because  they  both  began  with  the  results  of  a  quasi- 
paxeillel  local  analysis.  The  local  quasi-parallel  analyses,  as  discussed  in  Section  4.3 
yield  inconsistent  results.  Given  these  different  transients,  the  agreement  between 
the  two  PSE  results  is  some  indication  that  the  PSE  approximation  is  valid  for  this 
disturbance.  In  other  words,  the  viscous  terms  in  the  stability  equations  which  involve 
the  second  derivatives  of  the  shape  function  with  respect  to  and/or  axe  indeed 
negligible. 

Finally,  the  growth  rate  predicted  by  the  PSE  analyses  which  included  the  and 

derivatives  of  the  basic  state  axe  compared  directly  in  Figure  4.14  to  the  local 
analyses  which  neglected  the  same  terms.  The  large  nonpaxallel  effects  are  cleaxly 
evident. 

4.4.2  Leeward  Meridian 

Finally,  some  linear  PSE  results  are  shown  for  the  30  kHz  disturbance  along  the 
leeward  ray  that  was  investigated  earlier  with  the  local  analyses.  The  nonparallel 
flow  effects  seem  even  larger  here  than  they  did  along  the  windward  ray.  These 
effects  axe  shown  in  Figure  4.15.  Figure  4.15  compeires  the  results  from  the  linear 
PSE  analysis  which  included  the  and  derivatives  of  the  basic  state  to  the  local 
linear  results  which  neglected  the  same  derivatives.  (The  linear  PSE  analysis  used 
the  physical  velocity  components).  The  linear  PSE  and  local  analyses  axe  clearly  very 
different. 


4.5  Summary 

Extensive  stability  results  including  calculations  along  lines  other  than  the  windward 
and  leewaxd  meridians  could  not  be  obtauned  within  the  contract  period.  Nevertheless, 
our  calculations  support  the  contention  of  Stetson,  et  ed  (1992)  that  the  wavelength 
and  period  of  the  second  mode  scale  approximately  with  the  boundaxy  layer  thickness 
for  the  same  edge  velocity.  At  angle  of  attack,  the  windwaxd  boundary  layer  is 
much  thinner  (compressed  greater)  and  the  leewaxd  boundaxy  layer  is  much  thicker 
(compressed  less).  Hence,  the  wavelength  of  the  most  unstable  second  mode  in  the 
windward  boundary  layer  is  much  shorter  and  the  frequency  much  higher  than  at 
zero  angle  of  attack.  Conversely,  the  wavelength  of  the  most  unstable  mode  in  the 
leeward  boundaxy  layer  is  much  longer  and  the  frequency  much  lower  than  at  zero 
angle  of  attack. 
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Figure  4.2:  AFWAL  PNS  Basic  State  Profiles  of  the  Flow  Over  a  Sharp,  7®  Half-Angle 
Cone  at  2®  Angle  of  Attack.  Windward  meridian,  R  =  746  (interpolated).  Bow  shock 
at  =  26.7.  a)  Density,  b)  contravariant  velocity  component,  c)  Temperature, 
d)  contravaricint  velocity  component  x  100. 


Figiire  4.3:  AFWAL  PNS  Basic  State  Profiles  of  the  Flow  Over  a  Sharp,  7®  Half- Angle 
Cone  at  2®  Angle  of  Attack.  First  derivatives  with  respect  to  windward  meridian, 
R  =  746  (interpolated).  Bow  shock  at  =  26.7.  a)  Density,  b)  contravaxiant 
velocity  component,  c)  Temperatme.  d)  contravaxiant  velocity  component  x  100. 
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Figure  4.4:  AFWAL  PNS  Basic  State  Profiles  of  the  Flow  Over  a  Sharp,  7°  Half- Angle 
Cone  at  2°  Angle  of  Attack.  Second  derivatives  with  respect  to  windward  meridiaui, 
R  =  746  (interpolated).  Bow  shock  at  =  26.7.  a)  Density,  b)  contravariant 
velocity  component,  c)  Temperature,  d)  contravarizmt  velocity  component  x  100. 
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Figure  4.5:  AFWAL  PNS  Basic  State  Profiles  of  the  Flow  Over  a  Sharp,  7®  Half- Angle 
Cone  at  2®  Angle  of  Attack.  Leeward  meridian,  R  =  712  (interpolated).  Bow  shock 
at  =  37.44.  a)  Density,  b)  contravariant  velocity  component,  c)  Temperature, 
d)  contravariant  velocity  component  x  100. 
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Figure  4.6:  AFWAL  PNS  Basic  State  Profiles  of  the  Flow  Over  a  Sharp,  7°  Half- Angle 
Cone  at  2®  Angle  of  Attack.  First  derivatives  with  respect  to  leeward  meridian, 
R  =  712  (interpolated).  Bow  shock  at  =  37.44.  a)  Density,  b)  contravariant 
velocity  component,  c)  Temperature,  d)  contravariant  velocity  component  x  100. 
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Figure  4.7:  AFWAL  PNS  Basic  State  Profiles  of  the  Flow  Over  a  Sharp,  7®  Half- Angle 
Cone  at  2®  Angle  of  Attack.  Second  derivatives  with  respect  to  leeward  meridian, 
jR  =  712  (interpolated).  Bow  shock  at  =  37.44.  a)  Density,  b)  contravaxiant 
velocity  component,  c)  Temperature,  d)  contravariant  velocity  component  x  100. 


Figure  4.8:  Comparison  of  Spatial  Growth  Rates  of  70  kHz  Disturbance  Along  Wind¬ 
ward  Meridian.  /?  =  0,  length  scale  =  x^fy/R^.  □  Physical  velocity  components 
(parallel),  O  Physical  mass  fluxes  (parallel),  A  Physical  velocity  components  (quasi¬ 
parallel),  V  Physical  mass  fluxes  (quasi-parallel). 


Figure  4.9:  Comparison  of  Spatial  Growth  Rates  of  70  kHz  Disturbance  Along  Wind¬ 
ward  Meridian.  Physical  components,  locally  parallel  flow,  /8  =  0,  length  scale  = 
x^fy/Re^,  fifth  order  interpolants  in  Solid  line:  First  order  interpolants  in  and 
Dashed  line;  Third  order  interpolants  in  and  Dot-dashed  line;  Fifth  order 
interpolants  in  and 
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Figure  4.10:  Comparison  of  Spatial  Growth  Rates  of  70  kHz  Disturbance  Along  Wind¬ 
ward  Meridian.  Physical  components,  quasi-paradlel  flow,  =  0,  length  scale  = 
x^/\/Regi,  fifth  order  interpolants  in  □)  First  order  interpolants  in  and  O) 
Third  order  interpolants  in  and  A)  Fifth  order  interpolants  in  and 


Figure  4.11:  Comparison  of  Spatial  Growth  Rates  of  70  kHz  Disturbance  Along  Wind¬ 
ward  Meridian.  Physical  components,  quasi-parallel  flow,  13  =  0,  length  scale  = 
x^jy/Rtx^,  fifth  order  interpolants  in  Solid  line:  First  order  interpolants  in  and 
Dashed  line:  Third  order  interpolants  in  and  Dot-dashed  line:  Fifth  order 
interpolants  in  amd 
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Figure  4.12:  Comparison  of  Spatial  Growth  Rates  of  30  kHz  Disturbance  Along  Lee¬ 
ward  Meridian.  ^3  =  0,  length  scale  =  x^jyJRtx^.  □  Physical  velocity  components 
(parallel),  O  Physical  mass  fluxes  (parallel),  A  Physical  velocity  components  (quasi- 
parallel),  V  Physical  mass  fluxes  (quasi-parallel). 


Figure  4.13:  Comparison  of  Spatial  Growth  Rates  of  70  kHz  Disturbance  Along  Wind¬ 
ward  Meridian.  ,5  =  0,  length  scale  =  x^fy/R^.  Dotted  lines:  locally  pairallel  basic 
state.  Solid  lines:  nonparallel  basic  state.  □)  Physical  velocity  components,  O) 
Physical  mass  fluxes. 


55 


Figure  4.14;  Compaxison  of  Spatial  Growth  Rates  of  70  kHz  Disturbance  Along  Wind¬ 
ward  Meridian.  =  0,  length  scale  =  Solid  lines:  local  analyses.  Dot- 

dashed  lines:  linear  PSE  analyses.  □)  Physical  velocity  components,  O)  Physical 
mass  fluxes. 


Figure  4.15;  Comparison  of  Spatial  Growth  Rates  of  30  kHz  Disturbance  Along  Lee¬ 
ward  Meridian,  /d  =  0,  length  scale  =  x^ly/E^.  □)  Physical  velocity  components, 
local,  O)  Physical  mass  fluxes,  local,  <d )  Physical  velocity  components,  PSE. 
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Chapter  5 


Hypersonic  Flow  over  a  Blunt 
Cone 


Stetson  et  al.  (1984)(STDS)  has  csurried  out  stability  experiments  on  flow  over  blunt 
cones  at  freestream  Mach  number  of  8.0.  The  half  angle  of  the  cone  is  7“.  It  is  40 
inches  long  and  9.823  inches  in  base  diameter.  The  effect  of  nose  blimtness  on  stability 
was  studied  by  using  different  nose  radii  of  0.15,  0.25,  0.50,  and  0.70  inches. 
The  stability  of  the  mean  flow  in  one  of  the  experiments,  the  flow  over  the  blunt  cone 
with  a  0.15-inch  nose  tip,  is  analyzed  in  this  section.  In  this  experiment,  freestream 
stagnation  temperature  is  1350  deg  J?,  freestream  temperature  is  96  deg  /2(54.3  deg  if), 
freestream  velocity  (cailculated  from  freestream  Mach  number)  is  1181.67  m/s  and  the 
unit  Reynolds  number  is  2.5  x  10®//t.  The  physical  interfaces  and  insert  files  for  this 
problem  are  in  the  subdirectory  w.wpafb. 


5.1  Mean  Flow  Calculation 

The  basic  flow  has  been  calculated  using  two  different  codes:  a  Thin  Layer  Navier- 
Stokes  (TLNS)  code  developed  by  Esfahanian  (1991);  and  the  AFWAL  Parabolized 
Navier-Stokes  (PNS)  code.^ 

Esfahanizm’s  code  solves  the  Thin  Layer  Navier  Stokes  (TLNS)  equations  using 
the  second  order  accurate  central  difference  algorithm  of  Beam  and  Warming  (1978). 
It  cdso  uses  a  scheme  to  fit  the  bow  shock  as  a  discontinuity  in  the  flow  so  that  the 
rest  of  the  shock  layer  can  be  computed  accurately  with  no  numerical  oscillations  and 
minimal  artifici2d  viscosity.  Esfahanian  (1991)  h2ts  established  grid  convergence  of  the 
solution,  and  the  results  from  his  finest  grid  calculation  (1300  by  200)  are  used  here 
for  the  stability  analysis.  Some  representative  profiles  are  presented  below. 

The  AFWAL  PNS  code  is  briefly  discussed  in  Section  4.1.  This  code  uses  a  bow- 
shock  fitting  scheme  also,  and  has  the  option  of  using  the  Beam- Warming  algorithm 
or  the  Roe  flux-differenced  upwind  s  'heme  (Rajendran,  1989)  to  solve  the  equations. 
Different  runs  were  tried  using  each  algorithm  on  different  grids.  One  grid  had  251  and 
the  other  501  grid  points  in  the  radial  direction.  The  grid  with  the  greater  number 

^see  Section  4.1  for  references. 
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of  radial  grid  points  required  smaller  and,  hence,  more  mairching  steps  because  of 
the  practical  limit  on  the  CFL  number.  At  the  saime  time,  the  PNS  code  implements 
some  logic  to  enforce  a  minimum  step  size  so  that  departure  (exponenti<dly  divergent) 
solutions  eire  avoided.  This  minimum  step  size  requirement  makes  it  impossible  to 
arbitrarily  refine  the  grid  in  the  radial  direction  and,  hence,  more  difficult  to  compute 
accurate  solutions.  Finally,  only  runs  with  the  Beam- Warming  algorithm  completed 
successfully,  and  only  with  significant  levels  of  implicit  smoothing,  explicit  artificial 
viscosity,  and  implicit  matrix  conditioning.  The  runs  with  the  Roe  upwind  scheme 
did  not  converge. 

Figures  5.1,  5.2,  and  5.3  compare  the  TLNS  and  PNS  basic  state  profiles  and 
their  first  two  derivatives  with  respect  to  Here,  has  the  meaining  assigned  to 
it  in  each  code  respective  coordinate  system.  See  Figures  5.7  and  4.1.  While  the 
meaning  of  is  different  in  each  case  -  the  TLNS  body-intrinsic  coordinate  being 
normal  to  the  body  surface  and  the  AFWAL  coordinate  being  normal  to  the  body 
axis  -  the  error  in  comparing  them  this  way  is  small  because  the  cone  angle  is  only 
7°  and  the  boundary  layer  is  thin.  Since  cos(7®)  =  0.9925  and  sin(7°)  =  0.1219, 
one  has  —  0.9925^pjvs-  This  should  have  only  a  very  minor  effect  on  the 

comparison  of  the  scal<ir  quantities  such  as  the  density  eind  temperature.  However,  the 
velocity  components  are  also  different.  The  velocity  components  shown  for  the  TLNS 
solution  are  the  physical  components  tangent  and  normal  to  the  body  surface,  whereas 
those  of  the  PNS  solution  are  the  and  contravariant  (axial  an,’  ladial)  velocity 
components.  This  difference  also  has  some  very  minor  effect  on  the  comparison  which 
the  reader  should  be  aware  of. 

Some  errors  in  the  AFWAL  PNS  solutions  are  clearly  evident  in  the  derivatives 
of  the  2LxiaJ  velocity  component  neau'  the  wall  (Figures  5.2  b  and  5.3  b).  These  errors 
eire  most  likely  caused  by  the  high  levels  of  matrix  conditioning  required,  especiailly 
considering  what  the  conduioning  does  to  the  Parabolized  Navier-Stokes  equations. 
Stalnaker,  et  al.  (1986)  (p.  110)  show  that  the  matrix  conditioning  effectively  changes 
the  eigenvalues  of  the  inviscid  fl’  x  vector  Jacobians  by  a  smeill  amount  e  in  the 
subsonic  region  of  the  boundary  layer.  The  high  levels  of  artificial  viscosity  required 
may  also  contribute  to  the  errors  in  the  results. 

Some  of  the  differences  between  the  two  solutions  may  edso  be  attributable  to  the 
transport  properties  used  in  the  basic  state  calculations.  Although  the  same  trans¬ 
port  properties  were  used  in  the  stabAity  analyses,  we  do  not  know  what  transport 
properties  the  AFWAL  PNS  code  uses. 

As  v’’  .  oe  seen  later,  all  of  these  errors  and  differences  together  have  a  significant 
e'l^ect  on  the  stability  analysis.  This  is  especially  true  since  the  largest  errors  or 
aifferences  occur  near  the  edge  of  the  boundary  layer,  close  to  the  critical  layer  of  the 
disturbance. 

Finally,  some  additional  TLNS  results  are  snown  in  Figure  5.4.  Figure  5.4  com¬ 
pares  the  TLNS  predicted  temperature  along  the  cone  surface  with  measurements 
of  STDS  experiments.  Discrepancies  here  are  mainly  due  to  two  reasons.  First,  the 
adiabatic  wall  boundary  condition  is  used  in  the  calculations,  while  experimental  con¬ 
ditions  may  not  be  100  percent  adiabatic.  Secondly,  the  thermodynamic  model  used 
in  the  calculations  (e.g.,  Sutherland’s  Law  for  viscosity)  may  not  accurately  repre- 
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Figure  5.1:  Comparison  of  TLNS  and  AFWAL  PNS  Basic  State  Profiles  of  the  Flow 
Over  a  Blunt,  7°  Half-Angle  Cone.  s/Rjsi  =  175  (interpolated).  Bow  shock  at  « 
150.  Solid  line:  TLNS  solution.  Dotted  line:  AFWAL  PNS  solution,  a)  Density, 
b)  contravariant  velocity  component,  c)  Temperature,  d)  contravariant  velocity 
component  x  100. 
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a)  b) 


c)  d) 

Figure  5.2:  Compaiison  of  TLNS  2uid  AFWAL  PNS  Basic  State  Profiles  of  the  Flow 
Over  a  Blunt,  7®  Half-Angle  Cone.  s/Rff  =  175  (interpolated).  Bow  shock  at 
«  150.  First  derivatives  with  respect  to  Solid  line:  TLNS  solution.  Dot¬ 
ted  line:  AFWAL  PNS  solution,  a)  Density,  b)  contravariant  velocity  component, 
c)  Temperature,  d)  contravariant  velocity  component  x  100. 
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Figure  5.3:  Comparison  of  TLNS  and  AFWAL  PNS  Basic  State  Profiles  of  the  Flow 
Over  a  Blunt,  7°  Half- Angle  Cone.  s/R^  =  175  (interpolated).  Bow  shock  at  « 
150.  Second  derivatives  with  respect  to  Solid  line:  TLNS  solution.  Dotted 
line:  AFWAL  PNS  solution,  a)  Density,  b)  contravariant  velocity  component, 
c)  Temperature,  d)  contravariant  velocity  component  x  100. 
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S/R, 


Figure  5.4:  Temperature  Along  the  Cone  Surface. 

sent  the  experimental  properties.  Typiccd  mean-flow  velocity  (parallel  to  the  surface) 
profiles  are  compared  in  Figure  5.5.  Differences  near  the  wall  axe  due  to  inaccuracies 
in  experimental  measurements  near  the  wall  due  to  the  large  probe  size.  Figure  5.6 
shows  that  the  shock  shape  predicted  by  the  TLNS  solution  is  in  good  agreement 
with  the  experimental  result. 

Virtually  all  of  the  blunt  cone  stability  analyses  were  conducted  with  the  TLNS 
baisic  state  solution.  The  stability  analysis  of  the  AFWAL  PNS  blunt  cone  basic 
state  was  performed  mainly  to  test  and  validate  the  software  modules  developed  for 
the  studies  of  the  sharp  cone  at  angle  of  attack.  As  an  additional  benefit,  this  work 
also  provided  some  indication  of  the  relative  accuracy  of  the  eigenvalues  that  can  be 
expected  when  using  the  AFWAL  PNS  basic  state  profiles. 
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Figure  5.6:  Shock  Shape  Comparison  with  Billig’s  Correlation  and  the  STDS  Exper¬ 
iment. 
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5.2  Reference  Scales  for  Stability  Analysis 

Boundaxy-layer  edge  values  of  the  basic  state  vary  with  streamwise  distance  for  a 
blunt  cone,  unlike  those  for  a  sharp  cone.  Therefore,  instead  of  using  edge  values  (as 
in  the  experiments  of  Stetson  et  al.)  fixed  basic-state  freestreaun  values  of  temperature 
(T^  =  54.3  K)  and  pressure  (p^  =  .004082743  atm.)  axe  used  as  reference  values.  The 
station  at  s*/Rn  =  175,  where  Rn  is  the  radius  of  the  spherical  nose  (  s*  =  0.66675m). 
W21S  input  as  a  reference  station.  The  reference  length  sceile  rl,  is  the  boundaxy-layer 
length  scale 

=  2.8511  X  10-“  m 

at  this  reference  station,  s*//?jv  =  175.  Hence  the  squaxe-root  of  local  Reynolds 
number  R  =  is  related  to  the  length  scale  as  rl  =  s’^^/Rref  and  nondimensionaJ 
distance  s  =  s*/rl  is  related  to  R  as 

R  =  yJ^Rref 

The  Prandtl  number  is  specified  as  0.72,  the  same  value  used  in  basic  state  solution. 


5.3  Stability  Analysis  Coordinates 

Two  different  coordinate  systems  are  used  for  the  local  stability  «ind  PSE  analyses 
here.  They  are  axisymmetric  body-intrinsic  coordinates  =  s,  =  rj,  =  <f>). 
and  local  Caxtesiam  coordinates  aligned  with  the  cone  surface. 


Body-Intrinsic  Coordinates 

This  coordinate  system  (Figure  5.7)  is  used  when  icoord=l  is  specified  in  the  input 
file  Definitions.  Transformation  from  Cartesian  coordinates; 


=  s;  axe- length  along  the  body 
=  T}]  normal  disteince  from  the  body  surface 
azimuthal  angle 

is  given  by  the  following  equations. 


X  =  {xi,(s)  -  77sm(^c(s))} 
y  =  {r6(5)  +  j7Cos(^c(5))}  sin{<l>) 
z  =  {nis)  -I-  7/COs(^c(5))}  cos(^) 
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Figure  5.7:  Body- Intrinsic  Coordinate  System  Used  for  Stability  Analysis 

where  ri(5)  is  the  radius  of  the  body,  X6(s)  is  the  axial  distance,  and  tan{0c(5)}  is 
the  slope,  at  the  body  surface  {rf  =  0).  Since  the  basic  flow  is  axisymmetric,  and 
disturbances  axe  periodic  in  direction  with  a  wave  number  n^. 

Local  Cartesian  Coordinates 

This  coordinate  system  is  used  for  stability  amaJysis  when  icoord*0  is  specified.  Here, 
all  curvature  terms  and  cone  divergence  aure  neglected. 

X  =  s 

y  =v 

z  =  z 
dx^ 

w 

dx^ 

de 

dx^ 


ds 

dy 

drf 

dz 


66 


5.4  Dependent  Variables 

The  codes  LSH  and  PSH  provide  several  sets  of  disturbance  variables  to  be  used  in 
the  analysis.  The  results  presented  for  the  blunt  cone  axe  obtained  using  p,  f  and 
physical  contravaxiant  velocity  components  in  curvilinear  coordinates  as  dependent 
variables.  For  this  h5^ersonic  flow,  grid  resolution  is  required  not  only  near  the  wall, 
but  also  at  the  boimdaxy-layer  edge.  Results  here  axe  obtained  using  200  to  220  grid 
points  in  the  normal  direction  with  uniform  spacing.  A  mild  logarithmic  stretching 
is  used  in  some  cases. 


5.5  Linear  Stability  Analysis 

Different  families  of  eigensolutions  exist  for  hypersonic  boimdary  layers,  as  pointed 
out  by  Mack  (1986).  Hence  careful  guesses  of  eigenvalues  axe  required  for  local  eigen¬ 
value  calculations  to  converge  to  the  desired  (usually  unstable)  eigenvadue.  Three 
families  of  spatial  eigenvalues  for  two-dimensional  waves  at  a  fixed  streamwise  (^^) 
station  {s*/Rn  =  175,  il  =  2338.535)  axe  shown  in  Figure  5.8.  NonparaJlel  effects 
due  to  basic-state  streamwise  gradients  and  normal  velocity  axe  tadien  into  account  in 
these  calculations.  Spatial  growth  rates  a  =  —a,  and  wavespeed  aire  plotted  against 
wavenumber  q,  in  the  top  and  bottom  plots,  respectively.  The  family  Su  consists 
of  stable  first-mode  solutions  and  unstable  second-mode  solutions.  The  phase  speed 
of  this  family  approaches  1  —  as  the  wavenumber  approaches  zero.  Two  fami¬ 
lies  {Sd\,  Sd2)  are  damped  for  all  wavenumbers.  The  wavespeeds  of  these  families 
approach  1  -i-  as  the  wavenumber  decreases. 

Typical  amplitude  profiles  of  eigenfunctions  p,  T,  u  and  u,  for  an  unstable  second 
mode  axe  shown  in  Figure  5.9.  It  is  clear  that  large  gradients  exist  not  only  near 
the  wall  but  also  near  the  boundary  layer  edge  where  the  critical  layer  is.  Either 
homogeneous  conditions  or  asymptotic  conditions  for  disturbances  axe  used  at  the 
freestream  boundary.  Amplitude  profiles  axe  predicted  accurately  by  both  boundary 
conditions  but  small  oscillations  in  phase  near  rimax  are  observed  in  solutions  obtained 
with  homogeneous  boundau'y  conditions.  Note  that  the  phase  is  obtained  by  talcing 
the  inverse  tangent  of  the  ratio  between  real  and  imagineiry  parts  of  eigenfunctions. 
Errors  axe  not  visible  in  real  and  imaginary  parts  themselves  as  they  are  very  small 
near  the  outer  boundary.  But  errors  axe  amplified  enough  to  be  observed  when  the 
ratio  is  taken.  The  use  of  asymptotic  boundary  condition  gives  ph2ise  profiles  that 
are  free  of  these  oscillations  and  allows  the  use  of  a  shorter  domain  in  •q. 

Growth  rate  vs.  frequency  plots  for  disturbances  of  different  azimuthal  wavenum¬ 
bers  Ud,  at  the  reference  station  is  shown  in  Figure  5.10.  Basic-state  nonparaillel  terms 
are  included  in  the  calculations.  For  three-dimensional  (3D)  waves,  both  first  mode 
and  second  mode  disturbances  are  found  to  have  unstable  frequency  ranges.  For  the 
2D  wave,  the  first  mode  is  stable  (with  basic-state  nonpeirallel  terms  included).  For 
the  second-mode  instability,  2D  waves  (n^  =  0)  are  found  to  be  the  most  imstable. 

The  growth  rate  of  2D  second-mode  disturbances  computed  using  the  AFWAL 
PNS  code  is  compared  to  that  obtained  with  the  TLNS  solution  in  Figure  5.11.  Two 
different  curves  are  shown  for  the  AFWAL  PNS  solutions.  One  curve  corresponds  to 
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Figure  5.9:  Amplitude  Profiles  of  the  Eigensolution  at  Station  s*/Rn  =  175  and 
F  X  10«  =  82.7. 
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Figure  5.10:  Growth  Rate  Variation  with  Frequency  at  s*/Rn  =  175. 


the  AFWAL  PNS  solution  with  251  radial  grid  points,  and  the  other  corresponds  to 
the  AFWAL  PNS  solution  with  501  radial  grid  points.  All  of  the  stability  einalyses 
include  the  effects  of  the  transverse  curvature  of  the  cone  but  neglect  the  velocity 
component  (radial  or  normal)  and  the  derivatives  of  the  basic  state.  (The 
derivatives  of  the  basic  state  are,  of  course,  identically  zero  because  the  flow  is  axi- 
symmetric.)  Moreover,  the  analyses  use  homogeneous  boundciry  conditions  near  the 
shock. 

The  Reynolds  number  R  shown  on  the  abscissa  in  Figure  5.11  is  the  square  root 
of  the  Reynolds  number  based  on  the  free-stream  conditions  and  the  axe  length  s* 
measured  along  the  body  surface.*  Also,  the  growth  rate  in  Figme  5.11  has  been 
nondimensionalized  with  the  length  scale  rl  =  s^/Rref- 

The  first  thing  to  note  is  that  the  growth  rates  computed  with  the  different 
AFWAL  PNS  solutions  axe  more  stable  than  the  growth  rate  computed  with  the 
TLNS  solution.  Moreover,  the  errors  in  the  AFWAL  PNS  solutions  lead  to  quite 
significant  errors  in  the  growth  rates.  In  fact,  the  growth  rate  obtained  with  the  501 
point  AFWAL  PNS  solution  is  approximately  10%  low  at  the  local  maximum.  The 
analyses  of  the  AFWAL  PNS  solutions  also  piodict  that  the  loc^d  maximum  occurs 
at  a  station  slightly  closer  to  the  apex  of  the  cone.  Nevertheless,  the  growth  rate 
computed  from  the  AFWAL  PNS  solutions  appears  to  be  converging  to  the  result 
obtained  from  the  TLNS  solution. 

Most  of  the  differences  between  the  computed  growth  rates  axe,  thus,  likely  at¬ 
tributable  to  inaccuracies  in  the  basic  state  solutions.  However,  if  the  AFWAL  PNS 
grid  is  refined  even  more,  it  may  be  that  the  different  nonpaxallel  effects  in  the  slightly 
different  coordinate  systems  will  still  cause  the  growth  rates  to  be  slightly  different. 
This  issue  could  be  resolved  by  comparing  the  PSE  results.  In  the  interest  of  studying 
a  three-dimensional  basic  state  flow  (Section  4)  within  the  contract  period,  though, 
we  did  not  spend  any  additional  time  solving  the  PSE  for  the  AFWAL  PNS  basic 
state.  Instead,  the  TLNS  blunt  cone  solution  was  used  for  all  of  the  following  detailed 
local  and  PSE  analyses. 

Vaxiations  of  growth  rate  with  Reynolds  number  il  at  a  fixed  frequency  {F  x  10®  = 
82.7)  for  2D  second-mode  disturbances  axe  plotted  in  Figure  5.12  and  Figure  5.13. 
The  effect  of  including  cone  curvature  terms  is  to  decrease  the  growth  rate  from 
that  obtained  using  local  Cartesian  coordinates  (Figure  5.12).  Including  basic-state 
nonpaxallel  terms  shifts  the  growth  rate  curve  to  the  left,  resulting  in  a  lower  value  of 
R  at  branch  I  neutral  point  (Figure  5.13).  Both  types  of  boundary  condition  at  the 
fax-field,  asymptotic  and  homogeneous  conditions,  essenticdly  give  the  same  results. 

Next,  vaxiations  of  growth  rate  (Figure  5.14)  and  wave  speed  (Figure  5.15)  with 
Reynolds  number  R  for  different  values  of  frequency  parameter  F  are  presented.  Basic- 
state  nonpaxallel  terms  axe  included  in  these  calculations.  The  range  of  unstable 
frequencies  gets  lower  with  increaising  R  while  the  range  of  wave  speeds  increcises.  This 
is  in  ax:cord  with  the  fact  that  the  wavelength  of  unstable  modes  scales  with  boundeiry- 
layer  thickness  which  increases  with  Reynolds  number.  Corresponding  amplitude 

*The  results  from  the  PNS  coordinate  system  had  to  be  transformed  according  to  s/Rm  = 
arccos(l  —  x^/Rn)  for  0  <  x^/Rn  <  1  —  sin(^c)  and  s/Rs  =  f  -H  {x/Rn  —  (1  —  sin(^e)))  for 
x^/Rfi  >  1  —  sin(?e) 
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Figure  5.11:  Comparison  of  Computed  Spatial  Growth  Rates  as  a  Function  of  Surface 
Arc  Length  Using  Two  Different  Basic  State  Solutions.  F  =  82.7.  □)  TLNS  solution, 
200  wall  normal  points.  O)  AFWAL  PNS  solution,  501  radial  points.  A)  AFWAL 
PNS  solution,  251  radial  points. 


Figure  5.12:  Effect  of  Curvature  on  Second-Mode  Growth  Rate  Variation  with  R  at 
F  X  10®  =  82.7. 
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Figure  5.13;  Effect  of  Nonpaxallel  Basic  State  on  Second-Mode  Growth  Ra.te  Vauriation 
with  R  a.t  F  X  10®  =  82.7. 


ratio  curves  (N  factor  curves)  are  given  in  Figure  5.16.  The  maximum  N-factor 
according  to  this  linear  stability  theory  in  about  5.6  near  the  end  of  the  cone  (s*/ Ri^f  = 
235). 

Figure  5.17  gives  growth  rate  vs.  R  curves  for  3D  second-mode  waves  at  F=82.7 
and  for  different  azimuthal  wave  numbers  n^.  Difference  in  growth  rates  between 
different  oblique  modes  becomes  less  at  large  R.  This  is  due  to  the  fact  that  for  the 
fixed  the  wave  angle  approaches  zero  as  the  radius  of  the  cone  decreases.  Neutral 
stability  curves  for  these  2D  second-mode  disturbemces  including  nonparadlel  effects 
of  basic  state  axe  nlotted  in  Figure  5.18.  Minimum  critical  Reynolds  number  for  the 
second-mode  instability  is  predicted  to  be  ar  Rcru  =  1480.  This  value  is  lower  than 
that  observed  by  Stetson  et  al.  Note  that  the  Reynolds  number  in  their  experiments 
is  defined  in  terms  of  the  boundary-layer  edge  reference  value.  Reynolds  number  in 
our  case  can  be  related  to  theirs  by 

R  ~  Rstetton^J iJ^ooPoot^^  I  {UtPtiioa) 


5.6  Results  of  Linear  PSE 

Local  solutions  obtained  from  LSH,  taking  nonpar2dlel  effect  of  me<in  flow  into  account 
when  approriate,  axe  specified  as  initial  conditions  for  the  PSE  calculating  PSH.  In 
addition  to  streamwise  gradient  of  the  mean  flow,  nonpzurallel  effects  due  to  streamwise 
vau’iations  of  shape  function  /  and  amplitude  modulation  axe  properly  t2dcen  into 
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Figure  5.14:  Second-Mode  Growth  Rate  Vciriation  with  R  at  Different  Frequencies. 


Figure  5.15:  Wave  Speed  Vau'iation  with  R  at  Different  Frequencies. 


Figure  5.16:  N-Factor  Variation  with  R  at  Different  Frequencies. 


R 

Figure  5.17:  Second-Mode  Growth  Rate  Variation  with  R  for  Different  Azimuthal 
Wavenumbers. 
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Figure  5.18:  Neutral  Stability  Curve  for  Second-Mode  2D  Disturbances. 

account  in  PSE  analyses.  Effective  growth  rate  for  a  disturbance  variable  q'  is  defined 
as 

=  i2c{^^}-  Im{ai} 

where  Qj  is  the  imaginary  part  of  the  streamwise  wavenumber  and  s  is  the  body 
intrinsic  direction.  Due  to  nonparallel  effects  (of  basic-state  and  disturbance),  growth 
rates  of  different  variables  will  be  different.  These  will  depend  not  only  on  s  (R),  but 
also  on  normal  distance  form  the  wall  rj. 

Variation  of  growth  rates  based  on  maocimum  streamwise  disturbemce  velocity 
Umax  at  different  frequency  parameters  F,  is  presented  in  Figure  5.19.  Corresponding 
amplitude  ratios  (N  factors)  are  given  in  Figure  5.20.  The  difference  between  this 
PSE  result  and  that  of  LSH  (Figure  5.16)  is  larger  at  lower  Reynolds  numbers  R  (or 
higher  frequency  range)  where  nonparallel  effects  aure  the  larger.  For  each  frequency 
parameter  F,  the  location  of  branch  I  neutral  point  is  predicted  by  linear  PSE  at  a 
higher  Reynolds  number  than  that  by  local  analysis. 

In  experiments  of  Stetson  et  al.  (1984),  mass  flow  fluctuations  are  measured  for 
growth  rate.  In  our  calculations,  massflow  fluctuation  is  defined  as 

m  =  pu  =  pu-\-  pu 

where  u  is  the  disturbance  velocity  component  in  s-direction  of  the  stability  coordi¬ 
nates.  Figure  5.21  shows  amplitude  profile  and  massflow  fluctuation  and  its  growth 
rate  variation  with  normal  distance  q  at  a  fixed  station  R  =  2309  and  F  x  10®  =  82.7. 
In  Stetson  et  al.  (1984)  growth  rate  is  defined  in  nomenclature  as 

2  dA 

(J  - - 

Rstetson  dRstetson 
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Figure  5.19:  Vaxiation  of  Growth  Rate  of  Maximum  u  with  R  at  Different  Frequencies. 
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Figure  5.20:  Variation  of  Maximum  Amplitude  of  u  with  R  at  Different  Frequencies. 
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Figure  5.21:  Profiles  of  Amplitude  and  Growth  Rate  of  Massflow  Fluctuations  at 
R  =  2309,  F  X  10®  =  82.7. 
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where  Rstetson  is  based  on  edge  values  and  distance  measured  in  surface  direction 
s.  This  implies  that  differentiations  is  in  s  direction  holding  7/  constant.  On  the 
other  hand,  they  are  obtained  by  measuring  the  mciximum  in  maissflow  fluctuations 
which  occurs  near  the  boundary  layer  edge.  Therefore,  it  is  not  clear  whether  normal 
distance  from  wail  rj  in  our  notation  is  held  constant  in  their  growth  rate  evaluation. 
In  any  case,  since  the  growth  rate  varies  rapidly  in  rj  near  this  location,  slight  variation 
in  measurement  position  will  cause  large  variation  in  measured  growth  rate. 


Figure  5.22:  Variation  of  Growth  Rates  with  R  aX  F  x  10®  =  100. 


Variation  of  growth  rates  with  R,  based  on  maximum  massflow  fluctuations, 
(^max)  and  maximum  streamwise  velocity  (umox),  are  given  in  Figure  5.22.  Growth 
rate  obtained  from  LSH  is  also  shown  for  comparison.  Note  that  simply  measuring 
a  different  variable  i.e.,  m  instead  of  u,  results  in  a  shift  in  predicted  bremch  I  and 
branch  II  neutral  points  further  downstream  (higher  R).  Using  m  gives  higher  growth 
rates  than  using  u.  Amplitude  ratios  at  different  frequencies  are  plotted  vs.  R  in 
Figure  5.23.  The  maximum  N-factor  is  higher  than  that  predicted  with  Umax  and  the 
evolution  is  quite  different  high  values  of  F  (compare  with  Figure  5.20). 

Figure  5.24  shows  evolution  of  maximum  amplitude  for  different  variables.  Again, 
all  amplitudes  are  scaled  by  freestream  (preshock)  fixed  scales  as  opposed  to  exper¬ 
iments  where  boundaxy-layer  edge  values  are  used  as  a  reference.  Initial  maximum 
amplitude  for  u  aX  R  =  1520  is  specified  at  0.01  percent  of  mean  flow  freestreaim 
value.  This  gives  the  initial  amplitude  of  massflow  disturbance  of  about  0.07  %.  As 
observed  in  by  Stetson  et  al.,  velocity  fluctuations  are  much  smaller  than  temperature 
fluctuations  and  massflow  (or  density)  fluctuations.  At  the  laist  streamwise  station, 
the  maximum  amplitude  of  T  is  about  6.56  percent  <ind  that  of  m  is  about  1.55  per- 
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Figure  5.23:  Variation  of  Maximum  Amplitude  of  m  with  R  at  Different  Frequencies. 


cent.  The  slopes  of  the  curves  indicate  that  the  growth  rate  of  the  maximum  u  is 
slower  compared  to  that  of  the  other  three  variables. 


5.7  Results  of  Nonlinear  PSE 

Linear  PSE  calculations  indicate  that  the  use  of  different  variables  for  caJculating  (or 
measuring)  growth  rates  will  give  different  results.  It  also  shows  a  smedl  shift  in  branch 
I  neutral  point  downstream  due  to  nonparallel  effects.  But  growth  rates  obtained  by 
linear  PSE  are  still  higher  than  experimental  values,  and,  hence,  this  discrepancy 
caimot  be  accounted  for  by  nonparallel  effects  alone.  In  experiments,  the  harmonics 
of  the  dominant  frequency  are  observed  at  large  R,  indicating  nonlinear  effects  are 
present.  The  branch  I  is  also  further  downstream  of  that  predicted  by  either  local 
stability  analysis  or  linear  PSE  analysis.  In  this  section.  Nonlinear  PSE  simulations 
are  performed  to  explain  the  observance  of  higher  harmonics  in  experiments. 

5.7.1  2D  waves 

Nonlinear  evolution  of  a  2D  second-mode  disturbance  of  fundamental  frequency  F  x 
10®  =  82.7  is  first  investigated.  The  fimdamental  frequency  chosen  is  close  to  the 
dominant  (fundamental)  frequency  in  experimental  measurements.  The  initial  dis¬ 
turbance  at  F  X  10®  =  82.7  is  calculated  by  LSH  at  R=1809.4  (s=1400)  including 
basic-state  nonparallel  terms.  This  disturbance  is  specified  as  an  initial  mode  for 
PSH.  The  maximum  index  of  temporal  modes  is  modt  =  4,  i.e.,  first,  second,  and  third 
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Figure  5.24:  Evolution  of  Maximum  Amplitude  with  R  for  1 
F  X  10®  =  100. 


Figure  5.25:  Variation  of  Amplitude  of  u  with  R  for  the  Initial 
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harmonics  of  the  initial  mode  are  created  by  nonlinear  forcing  terms  and  followed  as 
integration  proceeds.  Amplitude  variation  of  the  fundamental  wave  and  its  harmonics 
is  given  in  Figure  5.25  to  Figure  5.27  for  three  different  initial  amplitudes.  Growth 
rate  variation  with  R  is  plotted  in  Figure  5.28  for  the  initial  u  amplitude  of  2  x  10"^. 


Figure  5.28:  Variation  of  Growth  Rates  of  m  with  R  for  the  Initial  u  Amplitude  of 
2  X  lO-'*. 


Linear  growth  of  the  second-mode  disturbance  for  the  saime  initial  amplitude  of 
10“^  is  also  given  for  comparison.  For  the  case  of  lowest  initial  amplitude,  nonlinear 
growth  rate  of  the  fundamental  wave  is  close  to  that  predicted  by  linear  PSE,  as 
expected.  The  main  feature  is  that  disturbances  tend  to  decay  at  some  value  of  R 
downstream  (Figure  5.29)  for  larger  initial  amplitudes.  This  decay  in  amplitude  is 
not  apparent  for  initial  amplitude  of  4“^.  Even  though  the  growth  rate  in  this  case  is 
decreasing  more  rapidly  than  other  case  (Figure  5.29),  large  amplitudes  of  different 
modes  cause  transition  to  begin  earlier  before  decay  or  saturation  can  be  observed. 
The  unstable  range  of  R  is  decreased  as  the  amplitude  is  increased,  branch  I  moving 
downstream  and  branch  II  moving  upstream. 

The  variation  of  skin  friction  coefficient  for  three-different  initial  amplitudes  is 
given  in  Figure  5.30.  The  rapid  rise  in  skin  friction  near  the  end  of  these  plots  can  be 
taken  as  the  beginning  of  the  transition  process.  This  shows  that  transition  location 
is  dependent  on  initial  conditions. 
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Figure  5.29:  Variation  of  Growth  Rate  of  m  for  Fundamental  Mode  with  R  for  Dif¬ 
ferent  Initial  u  Amplitude. 


5.7.2  3D  waves 

Since  the  actual  transition  process  involves  three  dimensional  structures,  nonlinear 
development  of  3D  modes  are  studied  here.  First,  nonline^  interaction  of  a  2D 
second-mode  disturbance  at  F  x  10*  =  82.7  and  a  3D  first-mode  disturbance  of  half 
that  frequency  is  investigated.  The  azimuthal  wavenumber  for  the  3D  disturb2uice  is 
specified  as  20.  This  value  is  kept  constant  in  the  integration  process,  in  accord  with 
the  irrotationality  condition.  The  initial  u  amplitudes  of  these  two  modes  are  set  at 
10“^  and  2  x  10”^.  No  parametric  excitation  of  the  subharmonic  wave  is  observed  for 
this  case.  In  fact,  the  growth  rate  (based  on  meissflow  fluctuation)  of  the  first  mode 
is  decreased  by  nonlinearity,  eventually  leading  to  its  decay  around  R=2580. 

In  the  above  calcidation,  higher  heirmonics  of  both  first  and  second  mode  distur¬ 
bances  axe  neglected,  except  for  the  distortion.  Moreover,  the  behavior  of  this  first 
mode  depends  on  its  obliqueness,  i.e.,  azimuthad  wavenumber.  Therefore,  a  sinailar 
run  is  performed,  with  the  wavenumber  increased  to  30,  and  keeping  the  next  order 
harmonics,  i.e.,  modes  (2,2),  (3,1),  (0,2),  and  (4,0)  in  the  interactions.  Evolution  of 
maximum  amplitude  of  u  fluctuations  for  this  case  is  given  in  Figure  5.33.  Here,  the 
subhaxmonic  mode  (1,1)  does  not  decay.  It  grows  but  the  growth  rate  is  lower  than 
that  of  the  second  mode  (2,0)  or  the  steady  vortex  (0,2)  or  oblique  higher  modes  (3,1) 
eind  (2,2).  The  fundamental  second  mode  amplitude  remrins  an  order  of  magnitude 
larger  than  those  of  the  rest  of  the  modes.  Hence,  unless  the  initial  amplitude  of  the 
3D  subharmonic  mode  is  significantly  larger  than  the  2D  wave,  rapid  rise  in  mean 
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Figure  5.30:  Variation  of  Mean  Skin  Friction  Coefficient  with  R  for  Different  Initial 
Amplitudes. 


skin  friction,  and  the  beginning  of  breakdown  is  due  to  distortion  by  the  2D  wave. 

Since  the  dominant  frequency  observed  in  experiments  is  not  as  low  as  the  first¬ 
mode  3D  wave  considered  above,  transition  process  occur  most  likely  due  to  3D 
second-mode  waves  and  their  harmonics.  Hence,  the  interaction  of  two  oblique  second¬ 
mode  disturbemces  of  the  same  initial  amplitude,  the  same  frequency  and  symmetric 
in  azimuthal  direction  is  considered  next.  The  azimuthal  wavenumber,  is  chosen 
to  be  16.  This  gives  a  wave  angle  of  about  40  degrees  at  the  initial  location  So  =  900. 
The  frequency  is  F  =  82.7  (u;  =  0.1934).  Maximum  index  for  temporal  modes  is  3  and 
maximum  index  for  a.zimuthal  wavenumber  is  4.  Therefore,  there  axe  10  modes  in  the 
calculations.  Amplitude  evolution  with  R  for  various  modes  are  shown  in  Figure  5.34. 
Three-dimensional  modes  are  found  to  grow  rapidly  near  the  end  of  the  cone,  filling 
up  the  spectrum.  In  particular,  steady  streamwise  vortices  (0,2)  emd  (0,4)  modes  are 
found  to  be  growing  very  rapidly. 

The  skin  friction  coefficient  for  the  mean  flow  is  given  in  the  next  figure.  Since 
steady  vortices  (0,2)  and  (0,4)  modes  contribute  to  streeimwise  component  of  the  wall 
shezu:  stress,  steady  skin  friction  is  periodic  in  azimuthal  direction.  The  maximum 
(peak  in  <l>  variation)  amplitudes  of  skin  friction  is  given  in  the  same  figure.  It  is 
evident  that  these  vortices  play  a  significant  role  in  the  rise  of  mean  skin  friction. 
The  meanflow  temperature  at  the  wall  is  set  to  zero  at  the  wall.  Near  this  location  of 
breakdown,  heat  transfer  normal  to  the  wall  is  found  to  increase  rapidly.  A  similar  run 
was  conducted  with  adiabatic  boundary  condition  for  steady  modes  at  the  wall,  giving 
similar  results.  In  this  case,  mean  wall  temperature,  instead  of  wall  heat  transfer, 
rises  rapidly  near  breakdown. 
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Figure  5.31:  Interaction  Between  2D  Second  Mode  and  a  3D  First  Mode  of  Half  Fre¬ 
quency;  Growth  Rates  of  Massflow  Fluctuations  Measured  at  the  Meiximum  Massflow 
Fluctuation  for  the  Fundamental  (2,0)  Mode. 
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Figure  5.32:  Interaction  Between  2D  Second  Mode  and  a  3D  First  Mode  of  Half  Fre¬ 
quency;  Amplitudes  of  Massflow  Fluctuations  Measured  at  the  Maximum  Massflow 
Fluctuation  for  the  Fundamental  (2,0)  Mode. 
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Figure  5.33:  Interaction  Between  2D  Second  Mode  and  a  3D  First  Mode  Including 
More  Harmonics;  Amplitudes  of  Massflow  Fluctuations  Measured  at  the  Maximum 
Mcissflow  Fluctuation  for  the  Fundamental  (2,0)  Mode. 
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Figure  5.34:  Nonlinear  Interaction  Between  Two  Oblique  Second-Mode  Waves:  Am¬ 
plitudes  of  Massflow  Fluctuations  Measured  at  the  Maximum  Massflow  Fluctuation 
for  the  Fundamental  (2,0)  Mode. 


Skin  Friction  Coofficiant 


Figure  5.35:  Nonlinear  Interaction  Between  Two  Oblique  Second-Mode  Waves:  Mean 
Skin  Friction  Coefficient  at  the  Wall. 
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Chapter  6 

Concluding  Remarks 


An  efficient  method  for  predicting  transition  in  high-speed  flows  over  fedrly  general 
bodies  has  been  developed.  The  major  objectives  set  out  in  the  contract  have  been 
achieved.  Instead  of  having  different  codes  for  different  problems,  we  have  developed 
two  codes  which  can  perform  two  functions,  the  local  linear  stability  analysis  (LSH) 
and  PSE  analysis  (PSH),  and  provide  for  physics  interf2w:es  to  different  problems. 
Interfaces  for  the  flow  over  a  flat  plate  and  hypersonic  flow  over  sharp  and  blunt  cones 
axe  built  in.  The  user’s  manual  describes  in  detail  how  to  set  up  new  problems.  The 
code  has  been  extensively  validated  and  used  to  study  the  stability  of  the  compressible 
flow  over  a  flat  plate  and  hypersonic  flows  over  sharp  amd  blunt  cones. 

Most  of  the  results  obtained  are  consistent  with  existing  theoretical  and  numerical 
results.  Typical  results  for  the  flat  plate  boundary  layer  agree  with  those  of  MjJik  et 
al.  (1990)  and  El-Hady  (1991).  Some  of  the  calculations  on  hypersonic  flows  over  the 
cone  explains  (partially)  some  observations  in  experiments.  The  results  revead  the 
sensitivity  of  growth  rates  measured  in  experiments  due  to  high  nonparallel  effects 
neair  the  critical  layer  where  the  growth  rate  is  typically  measured.  A  downstream 
shift  of  the  branch  I  neutral  curve  and  branch  II  neutral  points  due  to  nonparallel 
effects  is  also  apparent. 

The  appearance  of  hau:monics  of  the  fundamentad  second  mode  in  experiments  is 
simulated  by  both  2D  and  3D  nonlinear  PSE  calculations.  The  results  show  a  rapid 
Ailing  of  the  spectrum  near  breakdown  at  the  end  of  the  cone  and  saturation  of  the 
modes  due  to  nonlinearity.  These  charaicteristics  have  been  observed  in  sharp-cone 
experiments  of  Stetson  et  al(1983).  The  earlier  decay  of  disturbances  compared  to 
the  linear  results,  and,  hence,  the  shift  in  branch  II  neutral  point  upstream  (nairrow- 
ing  of  unstable  region)  due  to  nonlinearity  is  also  seen.  The  results  indicate  that 
transition  over  a  blunt  cone  is  most  likely  due  to  interaiction  of  oblique  second-mode 
disturbances.  Similar  nonlinear  phenomena  were  observed  in  results  on  a  sharp  cone 
by  Chang  and  Malik  (1993). 

The  results  of  applications  not  only  validate  both  method  and  code  but  adso  serve 
as  examples  for  the  user  in  solving  new  problems  as  described  in  the  user’s  manuad. 
In  aill  applications,  the  PSE  methodology  hats  been  found  to  be  a  useful  and  sound 
concept.  Input  models  for  nonlinear  calculations  still  require  more  experience  with 
receptivity  and  transition  in  high-speed  flows.  In  contrast  to  the  range  of  subsonic 
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speeds,  studies  on  stability  and  transition  at  high  speeds  are  currently  more  research 
than  routine.  While  the  development  and  utilization  of  the  full  capabilities  of  the 
PSE  methodology  will  require  some  time  and  effort,  the  combination  of  the  codes 
LSH  and  PSH  is  a  major  step  toward  an  efficient  tool  for  the  engineers,  researchers, 
and  scientists  to  analyze  amd  predict  the  transition  process  of  high-speed  flows  over 
fairly  general  bodies. 
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